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ABSTRACT 


In a variety of aeronautical applications, the flow around conical bodies at incidence is 
of interest. Such applications include, but are not limited to, highly maneuverable aircraft 
with delta wings, the aerospace plane and nose portions of spike inlets. 

For such conical bodies, starting at moderate angles of attack, the flow separates from 
the lee side, forming two vortices. Although the vortex lift contribution is highly desirable, 
as the angle of attack increases, the vortex system becomes asymmetric, and eventually 
the vortices breakdown. This causes problems with stability in all directions. Thus, some 
control of the separation process is necessary if the vortex lift is to be exploited at higher 
angles of attack. 

The theoretical model which is used in this analysis has three parts. First, the “single 
line- vortex” model is used within the framework of “slender body theory”, to compute the 
outer inviscid field for specified separation lines. Next, the three-dimensional boundary 
layer is represented by a momentum equation for the cross-flow, analogous to that for a 
plane boundary layer ; a von-Karman/Pohlhausen approximation is applied to solve this 
equation. The cross-flow separation for both laminar and turbulent layers is determined by 
matching the pressure at the upper and lower separation points. This iterative procedure 
yields a unique solution for the separation lines and consequently for the positions of the 
vortices and the vortex lift on the body. 

In the last part, control of separation is achieved by blowing tangentially from a slot 
located along a cone generator. It is found that for very small blowing coefficients, the 
separation can be postponed or suppressed completely (i.e., separation is moved all the way 
to the leewaxd generator), in which case the results from R.T. Jones’s theory are recovered. 
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NOMENCLATURE 


English letter symbols : 

a characteristic width of the body 

b characteristic thickness of the body 

bj half width of the wall jet 
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= \Ja 2 — b 2 , geometric parameter of the ellipse 

integration constant 

friction coefficient defined by eq(A5.11) 

wall friction coefficient for zero pressure gradient defined by eq( 10.20) 
total lift coefficient defined by eq(3.5) 
vortex lift coefficient 

pressure coefficient defined byeqs(3.2) and (10.21) 
blowing coefficient defined by eq(10.16) 
exponential function used in boundary layer analysis 
function defined by eq(2.3) 
metric (or Lame) coefficients defined in eqs(8.1) 
boundary layer shape factor 

7=1 

denotes the imaginary part of the complex quantity involved 

= r/27T, vortex strength 

(also used as a constant) 

spreading rate of the wall jet (constant) 

constant used in the wall jet analysis 

second pressure gradient parameter defined by eq(A5.14) 

characteristic dimension in the x - direction 

lift force 

Mach number 

exponent in the boundary layer growth eq(8.13) 

(also exponent in the wall jet analysis) 
normal force 
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static pressure 
dynamic pressure 

= (a 4- b)/ 2, local radius of the circular cone 

denotes the real part of the complex quantity involved 

Reynolds number 

cross-sectional area of the cone 

separation criterion (function) 

velocity component along the ior(- axes 

dimensionless velocity along the x or £ - axes 

velocity component along the y or 77 - axes 

dimensionless velocity along the 77 - axis 

velocity normal to the line vortex 

velocity component along z or £ - axes 

Cartesian coordinates fixed at the apex of the wing 

dimensionless Cartesian coordinates 

complex variable in the horizontal flat plate plane 


Greek letter symbols : 

a angle of attack 

7 vortex sheet strength 

T circulation of the line vortex 

8 wing semi-apex angle in the (x, z) plane 
(also boundary layer thickness) 

8 y displacement of the separation point 

from the leading edge along the y - axis 

8 Z displacement of the separation point 

from the leading edge along the z - axis 

81 displacement thickness in the direction of a cone generator 

82 displacement thickness in the circumferential direction 

e eddy viscosity 

e cone semi-apex angle in the (x, y) plane 
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complex variable in the vertical flat-plate plane 

(also coordinate in the direction normal to the surface) 

coordinate along the circumference of the cross-section 

complex variable in the circle plane 

[also angular coordinate (windward point taken as zero 

(also boundary layer momentum thickness) 

momentum thickness in the longitudinal direction 

momentum thickness in the circumferential direction 

momentum thicknesses due to the mutual effect 

of the longitudinal and circumferential flows 

dimensionless coordinate across the boundary layer 

first pressure gradient parameter defined by eq(A5.6) 

molecular viscosity 

kinematic viscosity 

coordinate along a cone generator 

3.14159... 

fluid density 

complex variable in the ellipse plane 
wall shear stress 

wall shear stress with zero pressure gradient 
disturbance velocity potential 
complex potential 

complex potential due to a line source distribution 
components of x* defined by eqs(6.3) and (6.4) 
disturbance stream function 


Subscripts : 


BL 

refers to the 

cf 

refers to the 

e 

refers to the 

j 

refers to the 

L 

refers to the 


boundary layer 
cross-flow plane 
non-viscous external flow 
wall jet 

laminar boundary layer 


m refers to the position in the jet profile where velocity is maximum 

s refers to the separation point 

T refers to the turbulent boundary layer 

tr refers to the transition point in the boundary layer 

oo refers to the undisturbed flow field 

0 refers to the solution for leading edge separation 

1 refers to the right vortex position 

Superscripts : 

a complex conjugate of a 

Abbreviations : 

BL boundary layer 

LE leading edge 

SPLS separation point on the lower surface 

SPUS separation point on the upper surface 


VI 


TABLE OF CONTENTS 


Abstract i 

Acknowledgement ii 

Nomenclature iii 

Table of Contents vii 

List of figures xi 

List of tables xiv 

1. PROLOGUE 1 

1.1 Motivation 1 

1.2 Objective 3 

1.3 Outline 4 

2. INTRODUCTION 6 

2.1 Previous theoretical work 6 

2.2 A comparison between the three models of vortex separation 7 

2.3 Previous experimental work 9 

2.4 Slender body theory 10 

2.5 Remarks on the assumption of “conical flow” 12 

2.6 Conformal mapping in the cross-plane 12 

3. A FLAT DELTA WING WITH ATTACHED FLOW 14 

3.1 The flow model 14 

3.2 Pressure distribution 15 

vii 



3.3 Lift 16 

4. A FLAT DELTA WING WITH LEADING-EDGE SEPARATION 17 

4.1 The flow model 17 

4.2 Boundary conditions 18 

4.3 Vortex position and strength 21 

4.4 Pressure distribution 21 

4.5 Lift 23 

5. A FLAT DELTA WING WITH DISPLACED SEPARATION LINES 26 

5.1 The flow model 26 

5.2 Boundary conditions 26 

5.3 Uniqueness of the solution 28 

5.4 Existence of the solution ' 28 

5.5 Vortex position and strength 29 

5.6 Pressure distribution 30 

5.7 Lift 31 

6. A DELTA WING WITH ELLIPTICAL CROSS-SECTION 32 

6.1 The flow model 32 

6.2 Boundary conditions 34 

6.3 Uniqueness of the solution 35 

6.4 Existence of the solution 36 

6.5 Vortex position and strength 37 

6.6 Pressure distribution 38 

viii 



6.7 Lift 38 

7. A CIRCULAR CONE 40 

7.1 The flow model 40 

7.2 Vortex position and strength 40 

7.3 Pressure distribution 41 

7.4 Lift 42 

7.5 Summary of inviscid results 42 

8. THE BOUNDARY LAYER ON A CIRCULAR CONE AT INCIDENCE ... 44 

8.1 Introduction 44 

8.2 The three-dimensional boundary layer 45 

8.3 The boundary layer equations, on a circular cone 47 

8.4 Laminar boundary layer 51 

8.5 Turbulent boundary layer 53 

9. DETERMINATION OF THE SEPARATION LINES ON A CIRCULAR CONE 55 

9.1 Viscous/Inviscid interaction 55 

9.2 Converged solutions for laminar and turbulent boundary layers 56 

9.3 Comparison with experiments 56 

9.4 Pressure distribution 58 

10. CONTROL OF SEPARATION BY BLOWING 59 

10.1 Introduction 59 

10.2 The flow model 60 

10.3 The wall jet equations 61 


IX 


10.4 Wall jet separation 64 

10.5 Converged solutions before and after blowing 66 

10.6 Pressure distribution 66 

10.7 Lift 67 

11. EPILOGUE 68 

11.1 Discussion 68 

11.2 Conclusions 70 

11.3 Recommendations for further research 71 

References 72 

Figures 77 

Tables 125 

Appendix 1 : Complex potential for an expanding ellipse 129 

Appendix 2 : Distance around the edge of an ellipse 132 

Appendix 3 : Evaluation of the derivative d\jda 134 

Appendix 4 : Three-dimensional boundary layer equations 137 

Appendix 5 : Solution of the boundary layer equation 142 

Appendix 6 : Program listings 157 


x 



LIST OF FIGURES 


PROLOGUE 

1.1 Model of F-5F at 40° angle of attack in Northrop water tunnel 78 

1.2 Vortex formation over a slender delta wing at incidence 79 

INTRODUCTION 

2.1 Three models representing vortex separation 80 

2.2 Conformal mapping in the cross-plane 81 

A FLAT DELTA WING WITH ATTACHED FLOW 

3.1 A flat delta wing at incidence 82 

3.2 Schematic of the cross-plane streamlines 82 

3.3 Pressure distribution for a/e = 1 83 

3.4 Lift versus relative incidence §4 

A FLAT DELTA WING WITH LEADING EDGE SEPARATION 

4.1 “Single line- vortex” model gg 

4.2 Locus of vortex positions gg 

4.3 Vortex strength versus relative incidence 87 

4.4 Pressure distribution for a/e = 1 gg 

4.5 Contours of integration for the normal force 89 

4.6 Lift versus relative incidence gg 

A FLAT DELTA WING WITH DISPLACED SEPARATION LINES 

5.1 “Single line- vortex” model gi 

5.2 Locus of vortex positions 92 

5.3 Vortex strength versus relative incidence 93 

5.4 Pressure distribution for a/e = 1 94 

5.5 Lift versus relative incidence gg 


xi 


5.6 Vortex strength and vortex lift vs. distance of separation from leading edge . . 96 
A DELTA WING WITH ELLIPTICAL CROSS-SECTION 

6.1 “Single line- vortex” model 97 

6.2 Schematic of the details in the cross-plane 98 

6.3 Locus of vortex positions for a 5% thick elliptical cross-section 99 

6.4 Vortex strength versus relative incidence 100 

6.5 Vortex solution boundaries 101 

6.6 Pressure distribution for ajs — \ 102 

6.7 Contours of integration for the normal force 103 

6.8 Lift versus relative incidence 104 

A CIRCULAR CONE 

7.1 “Single line- vortex” model 105 

7.2 Loci of vortex positions on a circular cone 106 

7.3 Pressure distribution for 9 a = 157° and a/s = 2 107 

7.4 Lift versus relative incidence 108 

THE BOUNDARY LAYER ON A CIRCULAR CONE AT INCIDENCE 

8.1 Secondary flow on a cone at angle of attack 109 

8.2 Streamline divergence producing a thinner boundary layer 109 

8.3 The mechanism of vortex coalescence within the boundary layer 110 

8.4 The coordinate system for the boundary layer analysis on the cone 110 

DETERMINATION OF THE SEPARATION LINES ON A CIRCULAR CONE . . . 

9.1 Schematic of boundary layer separation in the cross-plane Ill 

9.2 Flow chart for viscous/inviscid interaction 112 

9.3 Converged solutions for e = 5° and a = 30° 113 

9.4 Comparison of predicted separation with experiments 114 

9.5 Modified pressure distributions 115 

xii 



CONTROL OF SEPARATION BY BLOWING 

10.1 Schematic of controlled boundary layer separation 117 

10.2 Wall jet profile 118 

10.3 Converged solutions before and after blowing 119 

10.4 Blowing parameter versus angle of attack 121 

10.5 Pressure distributions before and after blowing 122 

10.6 Lift versus blowing 124 


xiii 


LIST OF TABLES 


6.1 Summary of scaling laws for displaced separation 126 

8.1 Analogy between two-dimensional and conical boundary layers 127 

10.1 Separation criteria for boundary layer and wall jet 128 


XIV 



1. PROLOGUE 


1.1 Motivation 

In a variety of aeronautical as well as aerospace applications, the flow around bodies 
of general conical shape at high angle of attack needs to be studied. Such applications 
include, but axe not limited to, highly maneuverable aircraft with delta wings such as the 
F-5F in fig(l.l), the aerospace plane and the nose portions of spike inlets. 

The class of shapes of current interest may be generalized as conical bodies of various 
cross-sections. For example, the nose portion of the fuselage of the F-5F can be approxi- 
mated by a circular cone, while the rear portion, being thinner and flatter, can be treated 
as a thin elliptical cone. The advantage of a semi-infinite cone considered here is the 
simple geometry and the assumption of conical flow, both of which simplify the analysis 
significantly. 

For such conical bodies, even at small to moderate angles of attack, the flow separates 
from the lee side. From this separation, fluid with high vorticity is convected upwards, 
away from the body surface, so that the resulting flow pattern is quite different from that 
of attached flow in which the vorticity is only appreciable in the boundary layer. 

In general, the flow pattern over a slender conical body goes through the following 
stages as the angle of attack increases : 

(i) At zero angle of attack the flow is axisymmetric. 

(ii) At very small angles of attack (0° < a < 6°) the flow is attached everywhere. 

(iii) At small angles of attack (6° < a < 20°) the flow first separates and a symmetric, 
steady pair of vortices exists. 
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(iv) At moderate angles of attack (20° < a < 45°) the symmetric vortex system yields to 
an asymmetric steady one (fig(l.l)), with two or more vortex cores. 

(v) At large angles of attack (45° < a < 70°) the asymmetric vortex system becomes 
unsteady and the vortices change locations randomly with time. 

(vi) As a -» 90° the vortices become highly mixed and form a turbulent wake (end portions 
of vortices in fig(l.l)). 

The present work is concerned with stage (iii) and only with the conical part of the 
vortices shown in fig(l.l). The numerical ranges for the angle of attack given in parentheses, 
although representative, are by no means absolute, since they have been determined from 
experiments with specific shapes, Reynolds numbers and other characteristics. 

Figure(1.2) shows the vortex formation over the leading edges of a slender delta wing. 
For highly swept-back configurations, usually a shear layer separates from each of two 
cone generators, one on each side, which roll up into a pair of vortices. These primary 
vortices cause additional lift to be generated on the lee side surface. The steep pressure 
gradient between the minimum of pressure and the primary separation line causes a new 
flow .separation, which usually takes the form of a small secondary vortex. The effect of 
these secondary vortices on the lift is usually small. 

The contribution of vortex lift in the low range of incidence is highly desirable. As 
the angle of attack increases, however, and the vortex system becomes first asymmetric, 
then unstable and uncontrollable, a large dependance on vortex lift may cause serious 
problems with directional, rolling and longitudinal stability. Therefore, if the formation of 
the vortices could somehow be enhanced or suppressed as necessary, controlled flow could 
be extended to higher angles of attack. 

The location and the strength of the vortices and, as a consequence, the vortex lift all 
depend on the location of the separation lines (as will be shown in later chapters). This 
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leads to the idea of controlling the location of separation as a means of controlling the 
vortex lift. 


1.2 Objective 


The purpose of the present work is three-fold : 

(i) First, to explore the influence of the position of separation on the vortex parameters 
(location, strength, lift). This is done through an inviscid analysis of the outer field for 
arbitrarily chosen separation lines. 

(ii) Second, to determine uniquely the separation line locations through a boundary layer 
(viscous) analysis. Thus, the ambiguity introduced in the first step is removed. 

(iii) Third, to control boundary layer separation by wall jet blowing. This also requires 
a viscous analysis and is based on the idea that a thin high-velocity layer of fluid ejected 
tangentially to the surface of the body reenergizes the boundary layer and makes it less 
susceptible to separation. 

As was mentioned in the previous section, at high angles of attack (typical of highly 
maneuverable aircraft), the problem is not to get lift (since a large component of the thrust 
produced by the engines is vertical), but rather, to get rid of any asymmetries present in 
the vortex system. Thus, some reduction in vortex lift as a result of blowing is justified. 

An alternate way to stabilize the vortices would be blowing from the apex along the 
axes of the vortices (brute force approach). However, as will be shown in later chapters, 
controlling the conditions which produce these vortices (i.e., boundary layer separation), 
is a more effective way to achieve our goal. This is indicated by the fact that very little 
tangential blowing produces very large changes in the vortex system. 
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1.3 Outline 


In chapter 2 a s ummar y of experimental and previous theoretical work is given, and 
a comparison is made of the various models currently in use for flows over conical bodies 
at incidence. The choice of the “single line-vortex” model as well as the use of “slender 
body theory” and the assumption of conical flow are also discussed. Finally, a conformal 
mapping sequence which allows simple transformations of various cross-section shapes is 
shown. 

In chapter 3 an account is given of the Jones model (ref.13) for a flat delta wing at 
incidence with attached flow. This is the most basic of all the models and provides the 
linear lift dependence on a. In chapter 4 the Brown and Michael solution (ref.19) for 
the separated flow past a flat delta wing is discussed. The separation is assumed to take 
place along the sharp leading edges. This is the simplest model from which the vortical 
(non-linear with a) contribution to the lift can be determined. In chapter 5 the influence 
on the flow geometry and the lift as the separation lines are moved inwards towards the 
leeward generator is determined. For a sharp leading edge this results in a singularity 
along the leading edges. In chapter 6 the singularity is removed by considering rounded 
leading edges (elliptical cross-section). In chapter 7 the inviscid analysis is extented to 
cones of circular cross-sections and compared with the results of Bryson (ref.20), but for 
various locations of the separation lines. 

In chapter 8 the three-dimensional boundary layer on the circular cone is solved by 
an extension of the Karman/Pohlhausen integral method to conical flow. In chapter 9 
an iterative viscous/inviscid interaction scheme is introduced which allows the prediction 
of the separation lines on the cone. Comparison with observed separation lines from 
experiments is also made. In chapter 10 tangential blowing into the boundary layer is 
introduced as a means of controlling the position of the separation lines and ultimately 
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the vortex lift. Finally in chapter 11 the limitations of the theory are discussed, 
conclusions are presented and suggestions for further research are given. 


the main 
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2. INTRODUCTION 


2.1 Previous Theoretical Work 

The fully attached flow past slender delta wings at incidence was first modelled by 
Jones (1946, ref.13), following earlier work with similar results by Munk (1924, refill) 
and Tsien (1938, ref.12). A little later Ward (1948, ref.14) completed the picture of the 
“slender body theory for attached flow” which was subsequently reviewed and extended 
by Adams and Sears (1953, ref.15). 

The three-dimensional separated flow past inclined bodies is currently represented by 
three well-established inviscid models which are described below (fig(2.1)). 

The first model is the “rolled-up core” established in 1957 by Mangier and Smith 
(ref.28). In this model, the inner turns of the rolled-up vortex sheet are represented by a 
single line-vortex, while a few turns on the outside of the spiral are represented explicitly 
in a numerical treatment (fig(2.1a)). 

The second model is the “multiple line-vortex” established in 1967 by Sacks, Lundberg 
and Hanson (ref.22) and it is derived in the following manner. If on the vortex sheet that 
springs from the leading edge of the wing, lines are drawn along which the circulation is 
constant (constant jump in the velocity potential A </?), these will also be streamlines of the 
mean flow. Each such line starts at a point on the leading edge and follows a helical path 
on the sheet, turning about the axis of the vortex as it proceeds downstream, thus dividing 
the sheet up into ribbons. The circulation about each ribbon is the same along the whole 
of its length and if it is allowed to condense into a line- vortex, a “multiple line- vortex” 
model is obtained (fig(2.1b)). 

The third one is the “single line-vortex” model. It is the simplest available model and 
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preceded those described above. It was finalized by Brown and Michael in 1954 (ref. 19) 
following earlier work by Legendre (ref.16), Adams (ref.17) and Edwards (ref.18). In this 
model, explicit representation of the outer turns of the spiral sheet is omitted, so that the 
cut, which in the first model connects the end of the vortex sheet with the concentrated 
vortex, now extends from the line-vortex to its associated separation line on the body 
(fig(2.1c)). A more detailed comparison of the three models is undertaken in the next 
section. 


2.2 A Comparison between the Three Models of Vortex Separation 

The arguments in this section follow those in refs.30 and 32 which should be consulted 
if a more detailed discussion is desired. 

The main advantage of the “rolled-up core” model of Mangier and Smith is, of course, 
its greater realism in describing separation. Thus, it is not surprising that it gives the 
closest approximation to the real flow in the infinite Reynolds number limit. Using the 
panel-method terminology, the “rolled-up core” model constitutes a higher order method 
than for ex am ple the “multiple line-vortex” model of Sacks et al, which means that it gives 
greater accuracy for a similar number of elements and has therefore the ability to predict 
a smooth behaviour of the flow. Thus, when many turns of a rolled-up configuration need 
to be represented, the “rolled-up core” model is the best one to use. 

Of course, there is a price for all these advantages, i.e., a greater programming effort. 
The computing time is also greater for the same number of elements than that required 
by the “multiple line- vortex” model. Also, the absence of any representation of secondary 
separation, may become important in some applications. 

The main advantage of the “multiple line- vortex” model is its flexibility, which makes it 
possible to use one program to calculate very different vortex structures with only minimal 
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changes. The approximation of the real flow in the infinite Reynolds number limit is also 
very good, and it is superior to the one given by the “single line-vortex” model but inferior 
to the one of the “rolled-up core” model. 

Unfortunately there are also complexities associated with the use of the “multiple line- 

vortex model. More computational storage space is required compared with the “rolled- 

up core” model, and a large number of elements is essential for accuracy. In addition, if 

many turns of a “rolled-up” configuration need to be represented, the calculation will be 

disrupted by vortices from adjacent turns pairing-off and rotating one around the other. 

Also, since the integration of ordinary differential equations, which is required in the 

streamwise direction in order to find the shapes of the line-vortices, may be an unstable 

process, the shapes sometimes become chaotic. When they do not become chaotic the 

shapes turn out to be helices with their pitch becoming smaller as the streamline gets 

closer to the axis of the vortex. It follows that a line-vortex starting near the apex of a 

delta wing should follow a helix of very small pitch, and such a helix requires very many 

elements to describe it with any realism. As more vortices are introduced to increase 

the accuracy of the solution, the closer to the apex the first one starts, thus making the 
problem worse. 

The main advantage of the “single line-vortex” model is its simplicity. This feature 
makes it especially attractive for initial investigations, and its use usually reveals the 
underlying structure of families of solutions of the more realistic models. Simplicity is also 
an important advantage when it becomes necessary to iterate the inviscid solution with a 
boundary layer solution in order to determine the separation lines. 

A disadvantage which arises with the use of the “single line-vortex” model is the 
inability to find solutions for very small relative incidences (a/e). For example, for the 
symmetric flow past a circular cone Bryson (ref.20) found no solutions with the vortex close 
to the separation line when (a/e) < 1.5 esc $. y whereas solutions have been found with the 
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“rolled-up core” model. In addition, since the vortex system is represented only globally 
in this model, the position of the vortices suffers in accuracy especially for asymmetric 
configurations. 

For the present work, the main purpose is to obtain a fast estimate of the velocity 
and pressure fields around the body which, when combined with a boundary layer analysis 
including the effects of blowing, will enable us to predict the separation lines. The “sin- 
gle line-vortex” model seems adequate for this purpose and will serve to demonstrate an 
approach that may be subsequently applied to the more elaborate models. 


2.3 Previous Experimental Work 

The experimental observations of separated flows on conical bodies, although limited, 
offer some very useful guidelines for the solutions which follow in the next chapters. 

Jorgensen (1957, ref.34) was the first to test cones with elliptical cross-section. He 
pointed out that there are distinct aerodynamic advantages to the use of elliptical cones, 
namely, that with their major axis horizontal, they develop greater lift and have higher 
lift-to-drag ratios than circular cones of the same fineness ratio and volume. However, his 
“lift coefficient versus angle of attack curves” are all linear, probably because the range of 
incidences tested was not high enough, so the vortex system either had not formed yet or 
was still too weak to affect the lift significantly. 

Rainbird, Crabbe and Jurewicz (1963, ref.35) experimented with circular cones in a 
water tunnel, while Schindel and Chamberlain (ref.39), and Friberg (refs.36,37) at M.I.T. 
tested circular and elliptic, two-dimensional and three-dimensional bodies. One of the 
interesting points of their results is the discovery of a secondary vortex system, similar 
to that shown in fig(1.2), above certain angles of attack. The observed positions of the 
separation lines from their experiments will be used for comparison with the predictions 
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of the present theory in chapter 9. 

Finally, Wood and Roberts (ref.64) at Stanford University showed that it is possible 
to control the vortex system by blowing tangentially from the leading edge of a slender 
wing, towards the leeward generator. Their work provides much of the motivation for the 
present analysis. 


2.4 Slender Body Theory 


From the mathematical point of view, “slender body theory” begins with the Prandtl- 
Glauert equation 


{ °° } dx 2 dy 2 dz 2 


( 2 . 1 ) 


which is valid for supersonic as well as subsonic Mach numbers. Equation(2.1) may be 
further approximated for slender elongated wings or bodies. It follows that, since the 
geometrical properties of the body or wing vary only slowly in the x - direction, the 
derivative d 2 <p/dx 2 must also be small. This argument can be made more rigorous (ref. 15) 
by introducing dimensionless coordinates 


X = lx * 

(2.2a) 

<< 

II 

o 

<Ci 

* 

(2.26) 

z = az* 

(2.2c) 
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where / is a characteristic length and a is a characteristic width of the body. If a function 
F is defined such that 


<?(x, y, z) = y m , z ’) 


(2.3) 


then eq(2.1) becomes 


(1 -ML) 


\ l 2 ) dx * 


d 2 F d 2 F d 2 F n 
+ dy-‘ + cl.- 1 _ ° 


(2.4) 


so that for sufficiently small values of the parameter (1 — M^)(a 2 // 2 ) the first term can be 
neglected. One thus obtains the Laplace equation for the cross-flow 


dV gV _ n 

dy 2 dz 2 


(2.5) 


However, the interpretation of slenderness is quite different for the various speed regimes. 
For Mo,, > 1 “slender” means that the wing lies well within the Mach cone from the 
apex. That was the reason why Ward (ref.14) limited his theory to pointed bodies and 
wings. Relatively blunt bodies and wings may qualify as “slender” at low supersonic 
speeds, whereas at hypersonic speeds eq(2.1) is not valid and the theory fails entirely for 
most practical shapes. For M*, < 1, on the other hand, the word “slender” becomes less 
restrictive, although we must keep in mind that eq(2.1) is not valid in the transonic regime. 
Since leading edge separation is essentially confined to highly swept wings, the application 
of “slender body theory” seems appropriate. 
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2.5 Remarks on the Assumption of “Conical Flow” 


The term “conical flow” implies that there is a point in the flow field, called the vertex 
of the flow, such that all physical quantities are constant along rays drawn from the vertex. 
The simplest example is the inviscid supersonic flow past a circular cone at zero incidence, 
with an attached shock. 

Strictly speaking, the flow over a conical body must be supersonic everywhere for the 
“conical flow” model to apply, since for subsonic flow the boundary conditions at infinity 
cannot be satisfied by a conical flow field. Nonetheless it has been observed that the 
subsonic flow past a slender conical body is approximately conical in a region downstream 
of the apex and well upstream of the trailing edge. This is due to the fact that at relative 
incidences (o/e) sufficient to cause separation, the circumferential pressure gradients are 
much larger than the axial pressure gradient caused by thickness and base effects. For a 
more thorough discussion on the subject of “conical flow” ref.7 should be consulted. 

The assumption of conical flow is facilitated in this model by two other assumptions. 
First, that the cone is of infinite length, thus avoiding the trailing edge region where the 
conicality assumption would break down. Second, the use of “slender body theory” which 
does not distinguish between subsonic and supersonic regimes since the first term in the 
Prandtl-Glauert equation is neglected, and there is no “upstream influence”. 


2.6 Conformal Mapping in the Cross-Plane 

Figure(2.2) shows the various relations which transform the cross-section of the cone 
from a flat plate (Brown and Michael solution), to an ellipse (Schindel solution) and finally 
to a circle (Bryson solution). The flat plate is most easily solved when transformed so that 
the vortex system is symmetrical with respect to a vertical plate. For the ellipse, the 
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easiest way is by a transformation to a circle and application of the circle theorem, which 
allows one immediately to write the complex potential in terms of the vortex system and 
its image. 
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3. A FLAT DELTA WING 
WITH ATTACHED FLOW 

3.1 The Flow Model 


This chapter presents the Jones solution for a flat delta wing, since this will be the 
departing as well as the destination point for the solutions which follow in the next chapters. 
It is the departing point since it excludes separation, and it is also the destination point of a 
separated configuration as the vortex system is being suppressed by shifting the separation 
line from the leading edge inwards, toward the leeward generator of the wing. Although 
only the flat delta wing is mentioned here, the lift coefficient based on the projected wing 
area is exactly the same regardless of cross-section, within the framework of “slender body 
theory without separation”. The configuration of the model is shown in fig(3.1). The 
flow pattern in any cross-plane is the familiar two-dimensional flow caused by a flat plate 
normal to a free stream with velocity u x a (fig(3.2)). However, the scale of the flow field 
increases continually along the x - axis, and this fact gives rise to the three-dimensionality 
of the problem. The potential function for the flow at the wing is given by 


V? = ±u <yo a\fa 2 - y 2 (3.1) 

where the positive sign is for the upper surface and the negative sign is for the lower 
surface. From eq(3.1) it may be seen that the gradient of the potential (i.e., the velocity), 
is singular at the leading edges (y = ±a). 
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3.2 Pressure Distribution 


The pressure coefficient is defined as 


C = (3.2) 

3oo 

where q^o is the free stream dynamic pressure. Following the analysis in ref. 13, the pressure 
distribution in a cross-section at a distance x from the apex may be expressed as a function 
of the relative incidence a/e 



whereas the difference in pressure between the upper and lower surfaces is 


AC? , g (<* 
e 2 y/a 2 — y 2 \e 


(3.4) 


The pressure distribution from eq(3.3) has been plotted in fig(3.3). The singularity at the 
leading edge is the result of an infinite suction there, as the flow tries to make a 180° turn 
from the lower to the upper surface. In a realistic description of the flow such a singularity 
cannot exist, and it is necessary to introduce a vortex sheet at the leading edge which 
feeds a vortex whose strength is such that the singularity is removed. In other words, a 
leading edge Kutta condition must be satisfied. This more realistic description of the flow 
was first provided by Brown and Michael (ref.19) and is shown in the following chapter. 
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3.3 Lift 


The lift coefficient is based on the projected area of the body (2 a ■ x/2), and is given 
by 



q 0 0 ax 


The Jones solution yields the classical result from “slender body theory” 


(3.5) 


% = 2tt- (3.6) 

e 2 e 

which shows that the lift grows linearly with angle of attack (fig(3.4)). The fact that the 
lift (i.e. the area between the two curves in fig(3.3)) is finite despite an infinite pressure 
peak at the leading edge should not be surprising, since the singularity is of the 1/^/y type, 
and becomes an infinite slope when integrated over y. 

Equation(3.6) verifies what is already known from experiments, i.e., that Ci/e 2 is a 
function of a/e, and suggests that a similar relationship should be sought for the more 
complex separated flow. Equation(3.6) will also be derived as a particular case from the 
more general flow configuration in chapter 4. 

Experiments show that the linear dependence of lift on the angle of attack is a fairly 
good approximation for small angles of incidence. As the angle of attack increases, however, 
the lift departs rapidly from the Jones value. In the next four chapters, an effort is made 
to capture this departure. It is done by acknowledging the fact that the flow separates at 
some angle of attack, and the Jones model is no longer a realistic representation of the 
flow field. 
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4 . 


A FLAT DELTA WING WITH 
LEADING EDGE SEPARATION 


4.1 The flow model 

This chapter discusses the Brown and Michael solution for the separated flow past a 
flat delta wing. The configuration of the model is shown in fig(4.1). The flow is assumed 
to separate along the two leading edges and to give rise to a pair of line-vortices. Since 
the strength of these vortices must grow in the x - direction for a conical flow field, they 
must be fed with vorticity from the leading edge. Otherwise Kelvin’s theorem would be 
violated. The connection between the leading edges and the vortices is achieved with plane 
vortex sheets emanating from the leading edges. This is the simplest possible way to model 
separation. 

A solution is now sought to satisfy eq(2.5) subject to the appropriate boundary con- 
ditions (section(4.2)). This is easily done by introducing the complex potential for the 
flow 


x =ip + i t/> (4.1) 

Note that although eq(2.5) is two-dimensional, the three-dimensional character of the 
problem will still enter through the boundary conditions. After using the transformation 
on the lower half of fig(2.2), the complex potential may be written as 


X(C) = -*UooaC-i*ln 


C~Ci 
C + Ci 


(4.2) 
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where k = T/2n is the vortex strength. Transforming eq(4.2) back to the physical plane 
gives 


x(*) 


—tu^aV z 2 — a 2 — ik In 


\fz* — a 2 — ^/z^ — a 2 
Vz 2 - a 2 + >/z? - a 2 


(4.3) 


where zi is the location of the right vortex, and z x is the complex conjugate of z x . The first 
term represents uniform flow past the plate (for small a) while the second term represents 
a vortex pair in the leeward side. 


4.2 Boundary Conditions 

The conditions that the solutions of eq(2.5) must satisfy are the following : 

(i) Tangency condition on the wing. This is automatically satisfied by choice of the 
complex potential. 

(ii) Separation condition on the wing. The separation line has to be specified since the 
present inviscid model is unable to predict it. In this chapter, the separation condition is 
simply the Kutta condition that the flow leave the plate tangentially at the leading edge 
and is most easily obtained in the £ - plane, where this condition requires the presence of 
a stagnation point at the origin. When transformed back to the physical plane it reads 


Upoa 

k 



(4.4) 


(iii) The disturbances must vanish at infinity. This condition is also satisfied automatically 
by the complex potential. 
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(iv) The fluid pressure must be continuous everywhere. This condition, however, cannot 
be met with the present model. The reason is that straight vortex sheets cannot be aligned 
with the flow. This difficulty could of course be circumvented by assuming curved vortex 
sheets, which would form part of a three-dimensional stream surface. The solution then 
would provide both the shape and the strength of the sheet (refs. 28— 33). However, the 
problem is greatly simplified by assuming straight feeding sheets and past experience has 
shown that such a model does capture the main features of the flow. The last condition 
needs therefore to be replaced by the following : 

(iv)' The vortex system (feeding sheet and concentrated vortex) must be force-free since 
only the wing and not the fluid can sustain forces. This requires that the force on the 
sheet be cancelled by an equal magnitude and opposite direction force on the vortex. The 
force on the vortex arises from its inclination to the local velocity vector, which in turn, 
derives partly from the free stream component along the x - axis and partly from the 
cross-flow velocity at its location. Thus, the force per unit length of the vortex may be 
written as 


-igu n T 


a 

-igu„ 7 - 


(4.5) 


where 


7 = 


= const 
ax 


(4.6) 


is the vortex sheet strength and 


u n = -Uoo£— + (v + iw)x 
a 


(4.7) 
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The force on the feeding sheet arises from the growth in T along the length of the vortex. 
The vector force per unit length of each filament representing the vorticity lying between 
x and x + dx is 


iQUool{.z\ ~ a) (4-8) 

Setting the vector sum of the two forces equal to zero according to the previous discussion 
and taking the complex conjugate of the resulting expression yields a condition for the 
induced velocity at the vortex location 


(v - tiw)i = ettoo — 1 j (4.9) 

This velocity may also be calculated by differentiating the complex potential in the physical 
plane (eq(4.3)) after the effect of the right vortex has been subtracted 


= Tz + 7=7, (410) 

Equating the right sides of eqs(4.9) and (4.10) yields the second relation between the two 
unknowns 


ik 


A — a 2 + yj{z\ - a 2 )(z\ - a 2 ) yj(z\ - a 2 )^ 2 - a 2 ) 


z\ 1 a 2 
z\ - a 2 + 2 zi(z 2 - a 2 ) 


= eu oo 



(4.11) 
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Equations(4.4) and (4.11) must be solved to determine the unknown 


quantities k and z x . 


o 


4.3 Vortex Position and Strength 


Solving the system of eqs(4.4) and (4.11) numerically by the Newton-Raphson tech- 
nique (refs.65,66) yields the position and the strength of the vortex as a function of the 
relative incidence. The vortex location is shown in fig(4.2). As may be seen, the vortex 
moves away from the surface and closer to the center-line of the wing for increasing a/e. 

Figure(4.3) shows that the vortex strength grows almost linearly with the relative 
incidence, the only departure from linearity occurring when the vortex system first appears, 
i.e., at very low angles of attack. 


4.4 Pressure Distribution 


The first-order expression for the pressure coefficient is given by (ref. 19) 



v 1 + xir 


ui 


(4.12) 


and its three-dimensional character is revealed by the inclusion of u which is dip/dx. The 
first term on the right is necessary when the coordinates are fixed on the wing and are 
tilted through the angle of attack, as is the case in the present analysis. Obviously, there 
is no contribution to the lift from this term, since it is exactly the same for both surfaces. 

The velocities may be computed by differentiating the complex potential 
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[ax ) [da ax J [da ) 


(4.13) 



(4.14) 


(4.15) 


where of course for this case of a flat wing w = 0 at the surface. 

Figure(4.4) shows the pressure distribution on the surface of the flat cross-section for 
a/e = 1.0, as compared with the corresponding pressure distribution in attached flow. 
The sing ular ity which appeared in the Jones solution has now been removed, since the 
flow no longer has to negotiate the sharp turn at the leading edge, but there is a pressure 
jump there due to the vortex sheet, equal to gu^T/x. This pressure jump is necessary to 
generate the force on the vortex sheet which balances the force on the vortex. 

The very low pressure region on the upper surface is the vortex signature, and its 
position corresponds approximately to the lateral location of the vortex. The peak suc- 
tion is an indication of the vortex strength, while the width of the suction is inversely 
proportional to the distance of the vortex from the surface. 

The difference between the solid lines and the dotted lines in fig(4.4) is, of course, the 
vortex lift. 
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4.5 Lift 


The normal force is most easily computed by calculating the change in downward 
momentum through an infinite plane perpendicular to the longitudinal axis x of the wing 
at the trailing edge (Trefftz plane). Thus 


N = 



dzdy 


(4.16) 


Note that dy/dz is the velocity component in a plane perpendicular to the wing surface, 
and therefore it contains the upwash contribution of the free stream. Integrating with 
respect to z produces a contour integral of the velocity potential 



(4.17) 


The contour is shown in fig(4.5) and includes the cuts connecting the separation points 
with the centers of the vortices. The vortices may be included in the body without affecting 
the normal force, since the forces on them cancel those on their feeding sheets. In terms 
of the complex potential 


N = -QUoo^t jy xdz + J xpdzj (4.18) 

Note that the z in the first integral is the complex variable in the physical plane while 
z in the second integral is the real variable in the direction normal to the wing surface. 
Since rp = 0 on the body and is single-valued on the vortices and the feeding sheets, the 
second integral vanishes. Furthermore, the function x( z ) is analytic in the field external 
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to the contour ; hence, the integral is independent of the path provided that it encloses 
the original contour. The simplest way to integrate eq(4.18) is by transforming it to the ( 
- plane (fig(4.5b)) 


The integral of the logarithm can be evaluated by deforming the contour into a large circle 
whose radius — ► oo. Since there are no singularities between the original contour and the 
large circle the integrals are equal. The remaining integration is done along the vertical 
line between the branch points and yields 


N = -QUcoft { / X(C) 


N = xga^u^a + gu^T^i + Ci) 

Transforming back to the z - plane gives 


(4.20) 


N = ■KQa 2 u\ 0 a + gu^T (\Jz\ - a 2 + \[z\ - a 2 ^ (4.21) 

or in dimensionless form 


z 


Equation(4.22) contains Cl because for small angles of attack the normal force can be 
taken equal to the lift. This is in agreement with the well known result that for a lightly 
loaded wing (small perturbation flow) the induced drag is a second-order quantity. The first 
term, being identical to the right side of eq(3.6), is the linear contribution from “slender 


% = 2 ^ + J Ll 

e 2 e auoo 


i - a 2 + \Jz\ - a 2 a 


(4.22) 
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body theory”, while the second term is the non-linear vortex lift. Figure(4.6) shows both 
components, as well as the total lift, as functions of the relative incidence. It may be seen, 
that the vortex lift is initially very small, but as the relative incidence increases it soon 
becomes the do mi nant term in the total lift. At a/e = 1.0 it has approximately the same 
magnitude as the Jones lift while at a/e = 2.0 it is approximately twice as large as the 
Jones lift. 

Brown and Michael have also carried out a second order approximation to eq(4.22). 
The analytical expression for this result is given by 

^ = 2tt| + 4.987 3 + 1.322 3 (4.23) 

This result is very similar to the expression derived by Smith (ref.29) using his “rolled-up 
core” model for thin slender wings in conical flow 


t= 2 *7 +4 - 9 (!r (4M) 

Thus, it was verified that Cl/z 2 is a function of the relative incidence a/e even for the 
case of vortical separation. 
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5 . 


A FLAT DELTA WING 

WITH DISPLACED SEPARATION 


5.1 The Flow Model 

The effect of vortex separation on the lift of a slender delta wing was examined in the 
previous chapter. The next step is to explore the effect of the location of the separation 
lines on the formation of the vortex system and consequently on the vortex lift. The easiest 
way to perform this task is through an inviscid analysis, in which the separation lines are 
selected arbitrarily . The flat delta wing offers once more the simplest geometry. 

The flow configuration is shown in fig(5.1) and is identical to the one used in the 
previous chapter, except that the separation lines have now been displaced a distance S y , 
still along generators but closer to the leeward generator of the wing. 

The complex potential is still given by eqs(4.2) and (4.3), respectively, for the trans- 
formed and physical plane. 


5.2 Boundary Conditions 


The requirement for separation from a point a — 6 y in the physical plane translates 

into 


dx 

*Jc. 


= 0 


(5.1) 


where 
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C = - 2a) 


(5.2) 


is the corresponding stagnation point in the transformed plane. The positive sign is used 
for separation from the upper surface while the negative sign is used for separation from the 
lower surface. Substituting eqs(4.2) and (5.2) into eq(5.1) and transforming the resulting 
expression back into the physical plane yields 


°° = + — == ( 5.31 

k — a 2 — t\JS y (6y - 2a) + ty /S v {6 y - 2a) 

where the signs in front of the square roots correspond to the flow separating from the 
upper surface (i.e., the positive sign is used in eq(5.2)). 

The procedure for balancing the forces on the vortex system is the same as that 
described in section (4.2). The force on the concentrated vortex is still expressed by 
eq(4.5) while the force on the vortex sheet is now 


-*0Uoo7(zi — a + 6 y ) (5.4) 

Setting again the vector s um of the two forces equal to zero and taking the complex 
conjugate gives an expression for the induced velocity at the location of the vortex 


( v - iw)i = e«oo 



Equating the right sides of eqs(4.10) and (5.5) now yields 


(5.5) 
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Z\ 

a 2 + yj{z\ - a 2 ){^ - a 2 ) 


Z\ 

- « 2 )(^i 2 - « 2 ) + ^ - 2a)(z 2 - a*) 


2 i l 1 a 2 

z\-a?- ty/6,(6 u - 2a) Wj{ - a*) 



(5.6) 


The numerical solution of eqs(5.3) and (5.6) gives k and z x in terms of a/e and 8 y /a. 


5.3 Uniqueness of the Solution 

An interesting result which occurred when the separation lines were forced away from 
the leading edges of the wing, was the appearance of two more families of solutions. They 
are discussed in some detail in section (6.3) ; here only the one which seems to agree with 
physical observations regarding the locus of the vortex positions for increasing relative 
incidence will be considered further. 


5.4 Existence of the Solution 

Another interesting aspect of the model with displaced separation lines is the difficulty 
in finding solutions for small a/e. The farther away the separation line is moved on the 
lower surface, the higher the minimum value of a/e for which vortex solutions first appear. 
This may be justified physically from observations of the actual flow over a flat delta wing. 
Since this flow separates at the leading edges, it is normal to expect some difficulty in the 
formation of the vortex system when the separation line is forced away from its natural 
position. Also, since in this model no account of the viscosity has been taken so far, it 
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must be concluded that it is the kinematics of the flow field which prevent the formation 
of the vortices at low angles of attack. 

The situation is quite different when the separation line is moved on the upper surface 
however. Then, there is no difficulty in finding solutions. The (a/ s ) Tnin remains zero as it 
was for the Brown and Michael solution. This may be explained by the fact, that although 
the separation line again shifts away from its natural position, it now moves towards the 
attached flow solution (i.e., the Jones solution which was discussed in chapter 3). 


5.5 Vortex Position and Strength 

Figure(5.2) shows the vortex location for various positions of the separation line. For 
a given separation line (i.e., constant 6 y /a ) the change in vortex location for increasing 
relative incidence resembles in general the Brown and Michael result, which corresponds 
to leading edge separation, except that it is displaced inward and toward the surface as the 
separation line moves inward. For a given angle of attack (i.e. constant a/e), on the other 
hand, the vortex is displaced toward the leeward generator of the wing. As the separation 
line approaches the center-line of the wing (i.e., ( 6 y /a ) — * 1), both curves collapse into 
the center of the cross-section, and the Jones solution is recovered. An interesting point 
illustrated in fig(5.2) is the sensitivity of the vortex position to very small displacements of 
the separation line. For a/e = 3, for example, it is seen that by displacing the separation 
only 1% causes a 10% shift of the y - vortex coordinate. 

From the numerical solution, the coordinates of the vortex were related to the Brown 
and Michael solution through the curve-fit approximate expressions given below 
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(5.8) 


The growth of the vortex strength with relative incidence is, as expected, similar to the 
Brown and Michael case (fig(5.3)), and it is very sensitive to small displacements of the 
separation line ; for a/e = 3, a 5% displacement of the separation reduces the vortex 
strength by approximately 25%. 

An approximate expression for the decay of vortex strength with distance of separation 
from the leading edge is given by (fig(5.6)) 


k ~ k o 



(5.9) 


where the Jones solution ( k = 0) is recovered for ( 6 y /a ) = 1. 


5.6 Pressure Distribution 

The analysis and formulation of section (4.4) are also valid in this case. Although 8 y 
does not show explicitly in the equations, it affects the solution for the vortex position 
and strength - as was discussed in the previous section - and as a result, the pressure and 
force on the wing. 

The pressure distribution for ( 6 y /a ) = 0.05 is shown in fig(5.4). It is very similar to 
the pressure distribution for the Brown and Michael solution but the singularity at the 
leading edge has been reintroduced. 
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5.7 Lift 


The vortex lift coefficient (fig(5.5)) is also related to that for the Brown and Michael 
solution by means of a similar approximate formula (fig(5.6)) 

Clv — Clv o (5.10) 

Both figs(5.5) and (5.6) clearly show the diminishing of the vortex lift contribution for 
increasing distance of the separation line from the leading edge, recovering the Jones lift 
in the limit as ( 6 y /a ) — ► 1. From a comparison of the exponents in eqs(5.9) and (5.10) as 
well as the two curves in fig(5.6), however, it may be seen that the vortex lift decays much 
faster than the strength of the vortex. This is so because the vortex lift is affected both 
by the strength and the position of the vortex. As the separation line moves inwards, the 
vortex gets closer to the wing surface so that the area that can benefit from the increased 
circulation also diminishes. In other words, the decay of the vortex lift is the result both 
of a diminishing vortex strength and of a diminishing area of lower pressure on the upper 
surface of the wing. 

One may argue that the practical value of this solution (i.e., a flat delta wing with 
displaced separation lines) is minimal because of the leading edge singularity ; however, 
it represents the limiting case for a delta wing of finite thickness with very small leading 
edge radius. 
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6 


A DELTA WING 

WITH ELLIPTICAL CROSS-SECTION 


6.1 The Flow Model 

There are two reasons that make rounded leading edges desirable. The first one is to 
get rid of the singularity which reappeared at the leading edge of the flat delta wing when 
the separation line was displaced ; a rounded leading edge eliminates the requirement of 
infinite acceleration and the resulting infinite velocity and pressure, although some suction 
will still exist. The second one is the necessity for a rounded edge in order to control 
separation by blowing. 

From the analytical point of view, the elliptical cross-section is the most convenient 
one to consider, since it may be easily related to both the flat plate as well as the circle 
by means of the Joukowski transformation (fig(2.2)). The flow configuration is shown in 
fig(6.1), the only new element from the previous one being the thickness. 

It is now easier to write the complex potential in the circle-plane (fig(2.2)), since the 
circle theorem allows one to write directly the contributions of the vortex system and its 
image 


Xcf(O) = -iu 0 0 a 



(9 - e x ) (ee l + r 2 )\ 
+ -&)) 


( 6 . 1 ) 


The first term on the right side represents uniform flow past the circle (the angle of attack 
has been assumed small), and the second term represents a vortex pair on the leeward side 
together with its image inside the circle. 


To this result we need to add the complex potential for an expanding ellipse of constant 
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axis ratio, in order to satisfy the tangency condition on the surface of the cone 


Xa = X*1 + Xi2 


( 6 . 2 ) 


where 


<t + \Ja 2 - c 2 
Xai = «oo be In 


(6.3) 


and 


Xj 2 = -UooeS {x ln2^/x(l - x) - 1 + 2 } 


(6.4) 


The derivation of x*i an d X «2 is given in appendix 1. Xa represents a source distribution 
on the horizontal plane of symmetry of the cone, directly related to the thickness. The 
total complex potential in the physical plane is of course 


X = Xc/ + X* ( 6 - 5 ) 

The velocity field is computed in the same manner as in section 4.4, except that now the 
evaluation of the derivative dx/da is more involved (see appendix 3). 
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6.2 Boundary Conditions 


The position of the separation point in the physical plane can be represented by 
cr 3 (a — S y , 8 Z ) (fig(6.2)). Under the Joukowski transformation a 3 goes into a point 9, in the 
9 - plane, given by 



- 8 y + %8 Z + y(a — 8y ) 2 + 2 8 z (a - 6 y )i - 6* - c 2 


Here 8 y and 8 Z are related by 


(6.6) 


^ + |=1 (6.7) 

since a 3 is a point of the elliptical cross-section. Requiring the presence of a stagnation 
point at <t 3 is equivalent to 


dXcf 



= 0 


( 6 . 8 ) 


or, from eq(6.1) 


U oo<* 

k 


9j + 29 A - R 2 R 2 + 29 A - 9\ 1 9] 

(d 3 + 9 1 )(9 a 9 l - R 2 ) (9 3 - 9 1 )(9 3 9 1 + R 2 ) J 9\ + R? 


(6.9) 


From the force balance, referring to fig(6.2), we may derive, in a similar manner to 
that described in chapters 4 and 5, an expression for the induced velocity at the vortex 
location 
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<7\ + \I<Tl— c2 01 + \/ <T l — C 


c 2 ) 2 + 4R 2 (a l + Ja 2 1 -c 2 )(a l + yJa 


1 

(2cr x , 8 y + i8 z ^ 

— £Uoo + 1 , 

1 ( 6 

<Ti + + yjcr 2 ~ c 2 + yjar\ - c 2 _ 

V a a ) 



The numerical solution of eqs(6.9) and (6.11) gives k and <J\ in terms of a/ £ and 9, 


6.3 Uniqueness of the solution 


Given the thickness of the wing and the position of the separation point, three solutions 
for the locus of the vortex positions for increasing a/e were found again (fig(6.3)) : 

In the first solution, the vortex moves farther form the wing surface and becomes 
stronger as the angle of attack increases. This is the only solution which agrees with 
experimental observations. 
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In the second solution, which appears in fig(6.3) as an extension of the first one, the 
vortex moves closer to the wing and its strength increases as the angle of attack increases. 

In the third solution, the vortex is under the wing and again moves farther away and 
becomes stronger as the angle of attack increases. 

The second solution disappears when the separation point is exactly at the leading 
edge. It should be noted that all three families of solutions exist also for the limiting cases 
of the flat delta wing (see section (5.2)), as well as for the circular cone. 

Although it might be interesting to investigate the question of stability for the second 
and third solutions, they are regarded as unrealistic and will not be considered further in 
the present analysis. 


6.4 Existence of the Solution 

As is shown in fig(6.3) there is a minimum value of the parameter a/e below which no 
solution exists. This is in agreement with experimental observations (refs.34-39), although 
the theoretical ( a/e) m i n may sometimes be larger than its corresponding value from ex- 
periments. The discrepancy between the experimental and theoretical values in this case 
results from the inability to satisfy the force balance for the vortex system due to the 
oversimplified representation of the vortex sheet. 

In fig(6.5), the vortex solution boundaries are shown as functions of thickness and 
separation location. One sees that (a/e) m ; n becomes smaller as the thickness of the wing 
diminishes. In the limiting case of a flat delta wing separation begins immediately for any 
a > 0 (provided that the separation is fixed at the leading edge). This is also in agreement 
with experimental observations (refs.34-39). 
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6.5 Vortex Position and Strength 


Approximate scaling laws similar to the ones shown in section (5.5), were also derived 
for an elliptical cross-section 10% thick 


2/1 - y io 


-fe )—(■-; 


( 6 . 12 ) 


z t - Zio 


(6.13) 


l -^ 1 ' 05 = *„(!- I ) 1 '” < 6 - W ) 

Here T) is the distance from the leading edge to the separation point along the surface, and 
the reference values are the ones corresponding to the flow separating from the leading edge. 
By comparison with those for the flat cross-section (see table(6.1)) it may be concluded 
that there is (almost) no variation in the vertical distance of the vortex from the wing 
surface as the wing acquires thickness, whereas the horizontal distance of the vortex from 
the center-line varies in the same manner as for the flat cross-section. As for the vortex 
strength, it appears to decay more slowly for the thick cross-section. Its growth with a/ £ 
is shown in fig(6.4) for a cross-section 5% thick and it is almost identical to the one for 
the flat cross-section with the same location of separation. 


k ~ k o 
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6.6 Pressure Distribution 


Equation(4.12) may once more be used to compute the pressure coefficient, except 
that w is no longer zero at the surface for the case of the elliptical cross-section. 

The pressure distribution for an elliptical cone 20% thick, is shown in fig(6.6). Note 
that the infinite suction singularity at the leading edge has been removed, but there is 
still some suction there. Comparing fig(6.6) with fig(5.4), one sees that the vortex suction 
is now larger (more negative peak). This is an indication of increasing vortex lift with 
thickness. 


6.7 Lift 

The normal force may be calculated in the same manner as for the flat delta wing. 
Equations(4.16) through (4.19) are still valid with the appropriate change of complex 
variable (z into a). The source term may be omitted from the complex potential since it is 
axisymmetric and therefore does not produce any downward momentum. The integration 
is performed again in the C - plane (fig(6.7)), as described in section (4.5), and the result 
in the present instance reads 


N = eUoJXCi + Ci) + 2 gu^ira ^ 


R 2 + -r - R - 


4 R, 


< R+ 4R / 


(6.15) 


Transforming back to the cr - plane 


N = 7r gu^a^a + 



a + b 
a — b 




(6.16) 
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or, in the usual dimensionless form 


Cl 


= 2x- + — 211— *((i + J<t\ - 7 + (1 - 

e u^ara a — bj y \ a — 


( 6 . 17 ) 


Figure(6.8) collects lift curves for all the configurations considered thus far. The 
discussion in section (5.7) regarding the various curves for different 6 y is valid for elliptical 
cross-sections as well. A comparison with the flat delta wing (see also table 6.1) clearly 
shows that the general trend is to achieve higher Cl for given a/e as thickness is added 
on the wing. 

The approximate scaling law showing the variation of vortex lift with 8 y is again 
similar to the one for the flat cross-section 


Clv — Clv o 



( 6 . 18 ) 


As was the case with the vortex strength, thickness causes the vortex lift to decay more 
slowly with shifting of the separation point toward the leeward generator. 
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7. A CIRCULAR CONE 


7.1 The Flow Model 

The circular cross-section may be looked at as one of the limiting cases of the elliptical 
cross-section, as b — * a, the other one being the flat plate (b — > 0, chapters 3, 4, 5). The flow 
configuration is shown in fig(7.1), and was first treated by Bryson (ref.20). The complex 
potential is given again by eq(6.5) except that the source term is now different 


X*i = UooRe In 9 (7.1) 

The difference is that the ellipse semi-minor axis has been replaced by the circle radius, and 
the plane source distribution has been replaced by a linearly growing source distribution 
along the axis of the cone. The logarithmic part could have been taken care of directly 
by the conformal transformation (see fig(2.2)). The transformation, however, would not 
change b into R. Thus, it is worth noting, that conformal mapping, although it comes very 
helpful in handling the cross-flow, cannot fully take into account three-dimensional effects 
such as the expansion of the body. 


7.2 Vortex Position and Strength 

The conditions and approach for the solution follow from the case of the elliptical 
cross-section with b = 1.0 and therefore they will not be repeated here. The only change 
that may be worth including in the computer programs is to replace the position of the 
separation point as given by (a — 6 y , 8 Z ) for the case of the ellipse, by a separation angle 
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Q a (the angle between the windward generator and the separation line). 

As was seen from fig(6.5) the minimum value of a/e for which a vortex solution first 
appears when the separation point is fixed at the leading edge (9, = 90°) is 6.22. Because 
this value is unrealistically high for practical applications, a separation angle close to the 
one observed in most experiments with circular cones (9, = 145°) was taken as a reference. 
For such angle, no simple formulae could be derived to relate the vortex coordinates, 
strength and lift for small excursions from this separation location as was done for the 
flat plate and the elliptical cross-section. It may be stated however, that again the vortex 
moves closer to the surface of the cone and becomes weaker as the separation point shifts 
towards the leeward generator. 

Figure(7.2) shows the domain of vortex solutions. The lower boundary of this domain 
is a function of both the location of separation and the relative incidence. It is worth noting 
that for small separation angles (9 3 ) vortex solutions cease to exist before the vortex reaches 
the surface of the cone. The upper boundary of the domain is the equivalent Foppl curve 
(i.e., the locus of the limiting vortex positions for high angles of attack, ref.20) for the case 
of a circular cone. 


7.3 Pressure Distribution 

The pressure distribution for a circular cone is shown in fig(7.3) for a/e = 2.0 this 
time, since the separation would have to be moved quite close to the leeward generator 
in order to get solution for a/e = 1.0 as in the previous cases. The same features may 
once more be identified (suction due to the vortex, and jump due to the vortex sheet). 
It is also worth noting that the Jones pressure distribution, shown in fig(7.3) with dotted 
line, is almost identical to the one for vortex separation, up to an angle of almost 100° 
from the windward point. The main difference, however, between the two cases is that the 
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Jones pressure distribution has only one adverse gradient while the pressure distribution 
for vortical separation has two adverse gradients since the two flows (starting respectively 
at the windward and leeward points) move toward each other. 


7.4 Lift 


From section (6.7) we may also 


get the lift coefficient for the limiting case of b = a 


£* £ UooCra £ 


(7.2) 


The above expression is plotted in fig(7.4). The lift curves do not go through zero as 
in the previous cases because of the absence of vortex solutions for small a/ s. Thus it is 
implied that at the point where they start they are connected with the Jones lift curve by 
a vertical straight line. This means that the vortex strength does not develop gradually 
from zero as in the flat cross-section case, but rather, it jumps into a certain starting value 
for the first a/e for which solutions are found. 


7.5 Summary of Inviscid Results 

t 

From the inviscid analysis in chapters 3 through 7 the following conclusions may be 
drawn : 

(i) The lift on conical bodies at incidence has two components ; the Jones lift and the 
vortex lift. The Jones lift is calculated assuming attached flow everywhere on the wing 
surface and grows linearly with angle of attack. The vortex lift is computed (in the present 
analysis) with the “single line-vortex” model and grows non-linearly with angle of attack. 
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The fact that the two lift components are decoupled suggests that blowing is a practical 
solution for changing the lift on the body without changing its attitude. 

(ii) As the separation lines are moved from the leading edges toward the center-line of the 
wing, vortex lift is suppressed and in the limit, as the separation lines coincide with the 
center-line, the Jones solution is recovered. This suggests that displacing the separation is 
indeed a viable mechanism for controlling vortex position and vortex lift. 

(iii) The vortex lift increases with increasing thickness of the wing (assuming always the 
same position of the separation lines). This advantage however cannot be realized at small 
angles of attack due to increased difficulty in finding solutions. 

So far, the separation lines have been chosen arbitrarily. In reality however, the 
position of separation must be determined through a viscous analysis. In other words, 
the velocity and pressure fields computed for the outer inviscid field, are introduced into 
the boundary layer equations ; integration of these equations yields two locations where 
the boundary layer leaves the surface, one oh each side of the hypothetical separation line 
which was arbitrarily chosen for the inviscid analysis. This procedure is undertaken in the 
next chapter for the circular cone. 
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8. THE BOUNDARY LAYER 

ON A CIRCULAR CONE AT INCIDENCE 


8.1 Introduction 

The main purpose of the viscous analysis is to predict the line(s) on the surface of the 
cone along which the boundary layer will separate. 

In general, the boundary layer on a cone goes through the following stages as the 
angle of attack increases (ref.52) : 

(i) At a = 0 it is similar to that on a semi-infinite flat plate or airfoil section and it 
may be studied by plane-flow methods. 

(ii) At small a it thickens at the top of the cone and thins at the bottom due to the 
circumferential flow induced by the angle of attack. 

(iii) At a/e ~ 0.5 an adverse circumferential pressure gradient first appears at the top 
of the cone. 

(iv) At some higher a, a separation bubble appears embedded at the base of the 
boundary layer growing in extent as a increases. 

(v) At a/e ~ 1.0 the boundary layer is no longer thin and the vortex bubble already 
existing at the top is in the process of coalescence into a symmetric pair of strong steady 
vortices. 

(vi) At a/e > 1 (i.e. a slender cone at very large incidence) the circumferential flow 
becomes similar to the plane flow about a cylinder and a von-Karman vortex street is shed 
at the top of the cone. 

In the present analysis we are concerned with stage (v). 
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8.2 The Three-Dimensioncil Boundary Layer 


Here are discussed briefly the properties that distinguish three-dimensional boundary 
layers from two-dimensional ones. 

(i) Secondary flow. 

In three-dimensional flow there are always pressure gradients at an angle to the main flow 
direction, providing a centrifugal force which distorts the outer flow streamlines. In the 
case of a cone at incidence for example, there is a circumferential pressure gradient while 
the main flow direction is almost longitudinal. Since the pressure is constant across a thin 
boundary layer, particles following a streamline within the layer are subject to the same 
circ umf erential pressure gradient as are those following the outer streamline. However, 
the bo un dary layer particles have lower inertia and tend to take a course conforming more 
closely to the direction of the circumferential pressure gradient as is shown in fig(8.1). 

(ii) Streamline divergence. 

The normal growth of a two-dimensional boundary layer is due to diffusion of vorticity 
(fig(8.2a)). In a three-dimensional boundary layer over a surface curved transversely to 
the direction of the flow (fig(8.2b)) it is necessary for the flow to spread itself over a 
progressively wider extent of surface as it grows. This spreading results in a thinner layer, 
than in the corresponding two-dimensional case. If the same velocity gradient is sustained 
between the surface and the outer flow in the two cases, the boundary layer on a cone will 
be thinner by \fy/Z than on a flat plate, resulting in a skin friction greater by y/Z (ref. 46), 
provided that equal lengths for the growth of the boundary layers are considered in the 
two cases. Stated differently, the cone boundary layer is similar to that on a flat plate with 
Reynolds number three times as great. 


- 45 - 


(iii) Separation. 

In plane flow separation occurs when a reverse-flow velocity profile appears, or equivalently 
when r = 0. In three-dimensional flow such a criterion fails to establish the separation lines 
because there is no way to decide which component of the shear stress is the important one 
to consider. However, it may be observed that at a separation line the wall stream surface 
bifurcates, and at a line of reattachment (if such occurs) the two stream surfaces join again 
at the wall. Thus, at the base of the boundary layer, there is embedded a distinct bubble 
that does not exchange fluid with the rest of the flow. It is therefore possible to generalize 
the definition of a separation region in three-dimensions as a bubble of fluid embedded in 
the boundary layer between the solid boundary and a stream surface meeting the body in 
a closed curve and containing a sheet pattern of vorticity. 

Of course, the kind of separation which is of interest for the present problem occurs 
when the embedded vortex sheet coalesces to form strong concentrated vortices. The 
mechanism of coalescence is described in ref.45. The vortex sheet is represented by a 
series of individual vortices as in fig(8.3), while the effect of the wall is represented by the 
image vortices below the wall. If the fluid above the wall imposes no additional constraint 
(i.e. boundary layer of locally infinite thickness), each vortex would move more to the left 
toward the separation point under the influence of the induced field of its image. Vortices 
initially near the separation point tend to remain fixed, however, because the layer is 
supposed to remain thin. Therefore, each vortex moving upstream tends to overtake the 
vortex ahead of it, and coalescence into a single strong vortex ensues. 
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8.3 The Boundary Layer Equations on a Circular Cone 


The three-dimensional boundary layer equations are given in appendix 4, written in 
a general system of orthogonal curvilinear coordinates. For a circular cone a coordinate 
system like the one described in section (A4.3) is obviously convenient. The geodesic 
coordinates are taken to be the cone generators while the geodesic parallels are the circles 
swept by the meridional angle, so the corresponding metric coefficient is the local radius 
of the cone (fig(8.4)). This gives 


^ = 1 (8-la) 

h 2 = R(i) = £sin£ ~ (8.16) 

h 3 = 1 (8.1c) 

The last equality in eq(8.1b) is validated by the assumption that the cone is slender. For 
the aforementioned coordinate system on a circular cone at an angle of attack, eqs(A4.3)— 
(A4.6) reduce to 

continuity 


s'du 1 dv dw 

“ + ? ar + 7a? + f ac =0 
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momentum in f - direction 
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momentum in 77 - direction 


dv v dv dv uv 

d£ e£ arj d( ( 


1 dp 1 dr„ 

dr) g 3C 


momentum in £ - direction 


(8.4) 



Furthermore, eqs(A4.7) and (A4.8) become 


(8.5) 
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Integrating eq(8.4) across the boimdary layer (i.e. from ( = 0 on the body surface to 


C — + 00 outside the boundary layer), one obtains 


/ ( dv v dv dv uv 
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1 dp\ 
get dr)) 



( 8 , 8 ) 


The normal velocity component w, can be replaced by 
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(8.9) 
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from the continuity eq(8.2). When this substitution is made and integration is carried out, 
eq(8.8) becomes 


rj u(v ‘ - v)d<: + 7^1 v( ' 1 ' ~ v)<K + h^ (v ‘ - v)di 

s 0 0 o 
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The displacement and momentum thicknesses 8 and 8 respectively are defined as follows 
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( 8 . 11 /) 
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Incorporating the above definitions into eq(8.10) gives the integral form of the boundary 
layer equation for the cross-flow 


9 d@22 dv - 

v e ~p— + (02 + 2022 ) + £U t V e (S\ + 20 2 2 ) 


+ £U e V e [{n + 2)021 - 2022 ] = 

Q 


(8.12) 


Here n is the exponent in the boundary layer growth expression 


S = k BL $ n 


(8.13) 


Equation(8.12) may also be written as 


9 00 22 ( dVm \ 

Ve ~d^ + [ Ve ~d^ + £UeVe ) ( S 2 + 29 ^ 

+ £ u e v e \ S \ ^2 4" {n -4- 2)02i — 2022] = ~ 1 " (8 14) 

e 

This equation is similar to the corresponding momentum equation for a two-dimensional 
boundary layer, the primary difference being the presence of the last term on the left side 
which contains the momentum thicknesses due to the interaction of the longitudinal and 
circumferential flows. 
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So far no assumptions have been made regarding the state of the boundary layer. 
Therefore, the above derivations are valid for both laminar and turbulent boundary layers, 
on condition that in the latter case u and v denote the time averages of the respective 
velocity components. The primary difference, however, between the two cases (i.e. laminar 
and turbulent) will be the rate of growth of the boundary layer (eq(8.13)). In the laminar 
case n = 0.5 while in the turbulent case n = 0.8 (refs. 40 and 41). 

The last term in eq(8.14) was evaluated numerically for several cases (a, e) and several 
locations (77). Its maximum contribution to the total value of the shear stress on the right 
side, occured when A,, = 0 (eq(A5.6)), and was approximately 13% for the laminar layer 
and 21% for the turbulent one. At separation (A„ = -12), its contribution was only 0.6% 
and 0.9% respectively for the two cases. Thus, it seems reasonable to neglect this term ; 
when this is done, eq(8.14) becomes 


^w +MSi+2 ^ ) (w +cu h I f (815) 

This result is exactly analogous to the corresponding equation for the two-dimensional 
boundary layer. Thus, the Karman/Pohlhausen method can be applied. The solution fol- 
lows immediatelly from the two-dimensional case, and the procedure is shown in appendix 
4. Table(8.1) illustrates the analogy between the various quantities involved in the two 
cases. 


8.4 Laminar Boundary Layer 


The separation criterion for the two-dimensional laminar boundary layer was estab- 
lished analytically by Pohlhausen and may be written in the equivalent conical terms as 
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Following the analysis in appendix 5, the above expression leads to eq(A5.32), which is 
repeated here 


'dV 


Vs 
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(8.17) 


U and V are the dimensionless external velocities 



Uoo 


(8.18a) 
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whereas El is defined by 


El = exp 



(8.19) 


For sufficiently slender bodies and small angles of attack, as has already been assumed in 
the inviscid solution, 


u e — =>• £7 — 1.0 


( 8 . 20 ) 
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so eqs(8.17) and (8.19) reduce to 


SC L = EZ'V- 6 



-0.334 


(8.17a) 


El = exp 


6e 



(8.19a) 


Thus, if the outer field is known, eqs(8.19a) and (8.17a) can be integrated. Starting 
from the reattachment point (77 = 180° for the usual range of angles of attack) the function 
SCl is computed until the point is reached for which SCi = —0.334. This will identify 
the upper separation point tj 3u . The integration is then carried out from the windward 
point (77 = 0°) until the separation criterion is satisfied again at some location. This will 
identify the lower separation point rj s i. The two points at which the boundary layer leaves 
the surface will, of course, enclose the point at which the vortex sheet emanates in the 
inviscid outer solution. 


8.5 Turbulent Boundary Layer 


The separation criterion for two-dimensional turbulent boundary layer was determined 
experimentally by Nikuradse, and may be written in conical terms as 




( 8 . 21 ) 


This result may be tranformed in a similar manner as for the laminar boundary layer 


(appendix 4) 
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where Ej is now defined by 


E? — exp 
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If one assumes as in the laminar case that U ~ 1, eqs(8.22) and (8.23) reduce to 
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The integration of the turbulent boundary layer equations is carried out in the same 
way as for the laminar boundary layer, assuming that the boundary layer is turbulent 
from its start. For a given cone (e) and angle of attack (a) the upper separation point is 
almost the same as for the laminar case, since the boundary layer which develops from the 
upper reattachment point has very little space to travel. The lower separation point, on 
the other hand, will occur at a larger distance rj sl since the turbulent boundary layer can 
progress farther into an adverse pressure gradient before separating, due to its increased 
momentum near the surface. 
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9 


DETERMINATION OF THE SEPARATION 
LINES ON THE CIRCULAR CONE 


9.1 Viscous/Inviscid Interaction 

The boundary layer analysis in the previous chapter was made possible by assuming 
that the term containing the momentum thicknesses due to the interaction of the longi- 
tudinal and circumferential flows was negligible compared with the rest of the terms in 
the momentum equation. Thus, for both the laminar and the turbulent case eq(8.14) was 
reduced to eq(8.15) which is similar to the momentum equation for the two-dimensional 
boundary layer (table(8.1)). 

On the other hand, the derivation of eq(8.14) in the first place, was made possible 
by the use of the simplified coordinate system described in section(8.3), which led to the 
simple expressions for the metric coefficients of eqs(8.1). 

Unfortunately for an elliptical cone only one of these coefficients is unity, while the 
other two, when expressed in terms of the local coordinates, contain hyperbolic and trigono- 
metric functions which lead to more complicated form for the boundary layer equations. 
For this reason, the viscous analysis is restricted to circular cones only. 

The boundary layer separation in the cross-plane is sketched in fig(9.1). The matching 
of the viscous and inviscid flow fields is illustrated in fig(9.2), and is described below : 

Velocity distributions as functions of the angle 0 around the circular cross-section of 
the cone calculated (for a = b = R = 1) using the “single line- vortex” model are introduced 
into eqs(8.17a) and (8.19a) for the laminar case and into eqs(8.22a) and (8.23a) for the 
turbulent case. The integrations are carried out numerically by the Romberg method 
(refs. 65, 66). First, the starting point is taken at the upper reattachment point (6 = 180° 
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if a is large enough) and proceeding clockwise (for the right side of the cone) the point 
where the top boundary layer leaves the surface is identified. Similarly, starting at 0 = 0° 
and proceeding counterclockwise (for the right side of the cone again), the point where 
the boundary layer leaves the surface is identified when SCi = —0.334 or SCt = —3.75 
depending on the state of the boundary layer. 

For a given cone geometry and angle of attack, the only acceptable solution (in terms 
of the assumed separation angle) is the one which yields the same pressures at both points 
where the boundary layer leaves the surface. This assumption is justified by the fact that 
the two streamlines which separate from the surface form a bubble (fig(9.1)) inside which 
there is no flow and therefore the pressure must be uniform. 


9.2 Converged Solutions for Laminar and Turbulent Boundary Layers 


Fig(9.3) shows the converged solutions for a cone with e = 5° at a = 30°, for laminar 
and turbulent boundary layers. It may be seen, that the main difference between the two 
cases is the location of the lower separation. As was expected, when the boundary layer 
is turbulent, separation is delayed until a larger angle. The locations of the upper and 
inviscid separations as well as the vortex positions are almost identical for the two cases. 


9.3 Comparison with Experiments 


Friberg (ref.37) performed several experiments with circular cones in which the exter- 
nal flow was subsonic and the boundary layer was laminar. He was able to fit his observed 
separation lines reasonably well by the formula 
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= 90 + (73 - V51.4a-450)(0.76 + 0.024s) 


(9.1) 


for 10 < a < 30 and 5 < e < 20 


Jorgensen (ref. 34), on the other hand, with his experiments in which the external 
flow was supersonic and the boundary layer turbulent, revealed a separation angle of 147°, 
which is also the angle that Bryson (ref. 20) uses in his model. 

In fig(9.4) the experimental results from both references are shown together with 
predictions of the present theory. The flat part which is common to all the curves in the 
low range of angles of attack represents attached flow (no vortex solutions exist in this 
range). At a ~ 5° which corresponds to a/e ~ 1 separation first takes place and all the 
separation angles change rapidly as a increases. Finally, at a ~ 15° which corresponds 

to a/e ~ 3, each separation angle reaches a limiting value which remains constant as a 
increases. 

The agreement of the theoretical predictions with experimentally determined points 

is excellent. Most points seem to fall on the lower separation curve for the turbulent 
boundary layer. 

Another feature shown by Friberg’s experiments is that, although there is clearly a 
trend of the separation lines to move windward as a increases, surprisingly, there are 
exceptions such as the last point (a = 30°) for the e = 5° cone in fig(9.4). 

Lastly, it must be noted that Friberg’s experiments showed a separation angle which 
is a function of both a and £ and not only of their ratio a/e, while the present theory 
shows that the separation angle is almost a unique function of a/e. 
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9.4 Pressure Distribution 


The modified pressure distributions for a/e = 2 are shown in figs(9.5a) and (9.5b) 
respectively for the laminar and turbulent boundary layers. The flat portion of these curves 
represents the separation bubble where the pressure is required to be uniform. Although 
the initial inviscid pressure distributions are almost identical since the inviscid separation 
is almost the same (0, = 157° and 9, = 159°) for the two cases, the modified curves which 
result from the inclusion of the boundary layer show two characteristic differences. First, 
the separation bubble is larger for the laminar case (see also fig(9.3)), and second, the 
vortex suction is more pronounced for the turbulent case. 
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10 


CONTROL OF SEPARATION 
BY BLOWING 


10.1 Introduction 

So far it has been shown that the boundary layer on a circular cone at incidence, as 
it develops from the windward stagnation line towards the leeward generator, will sooner 
(when it is laminar) or later (when it is turbulent) separate due to the adverse pressure 
gradient which encounters. It is possible, however, to postpone this separation, by replacing 
the natural boundary layer with a turbulent wall jet. The increased (due to the jet) 
momentum near the surface reenergizes the boundary layer, thus allowing the viscous flow 
to remain attached for a larger arc. 

The mechanism of delaying the boundary layer separation through blowing is sketched 
in fig(lO.l). The wadi jet changes the location of separation, and this in turn changes all 
the vortex parameters (position, strength and lift). In other words, blowing changes the 
entire (inviscid) outer flow field by changing the conditions which generate this field. 

Although the behavior of a wall jet flowing around a curved surface has been the 
subject of study for almost two centuries, the idea of using a thin, high velocity, tangential 
jet of fluid to control the location of separation on wings with rounded edges is relatively 
new. Wood and Roberts (ref.64) have recently examined the practicality of such a scheme 
by performing a wind tunnel experiment in which a wall jet was used, as a cross-flow plane 
device, to control the separation and hence the positions of the associated vortices on a 
conical delta wing. 

The analysis in this chapter follows after Roberts (ref.63). Although in our case there 
is an external flow, the jet velocity is assumed to be much higher than the velocity of 


-59 - 



the outer field (i.e. Vj > v e ). Therefore the jet will be treated as issuing into quiescent 
surroundings. In addition, since the thickness of the boundary layer and the width of the 
jet are small compared to the local radius of the cone (i.e., 6/R,bj/R <C 1) curvature 
effects will also be neglected. 


10.2 The Flow Model 

The profile of the wall jet is shown in fig(10.2). The jet consists of two parts ; an 
inner flow adjacent to the wall having a highly non-linear velocity profile characteristic of 
a turbulent wall flow, and an outer flow having a velocity profile typical of a free turbulent 
plane jet. The jet emerges from a point source into the fluid and spreads, increasing its 
width an d decreasing its velocity due to turbulent diffusion in the jet and friction at the 
surface. At a distance 77 downstream of the jet exit the velocity Vj can be expressed as 


Vj = v m (rj)f (10.1) 

where v m is the maximum velocity, occuring at £ = C m( 7 7) > an d bj = bj(ri) is the half width 
of the jet (at which point vj = v m /2). The velocity profile in the outer flow (( > Cm) is 
assumed to take the form 


vj ~ v m sech 2 


M ~ Cm) 

. (1 — Cm) . 


for C > Cm 


( 10 . 2 ) 


This velocity profile is suggested by the classical free jet solution by Tollmien, modified to 
give Vj = v m at C = Cm- The constant kj is determined such that Vj = v m /2 at C = bj. 
Thus 
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= 0.8814 


(10.3) 


kj = tan} i 1 f -^= 

Va/2 

The velocity profile for the inner flow is assumed to depend on the variable (C/Cm) 1/n 
(where n = 7), as suggested by turbulent wall flow, and is chosen to give a maximum value 
Vj = v m at C = Cm- Thus 


Vj = v m 



(10.4) 


The value of Cm is determined by matching the second derivative of the velocity profiles 
given by eqs(10.2) and (10.4). The result is written 


Cm = &i( 1 4- fcjn) -1 


(10.5) 


10.3 The Wall Jet Equations 

In addition to the approximations mentioned in section (10.1) (i.e. that v, > v e and 
bj/R < 1), the assumption is also made that the contribution of the shear stress at the 
wall and the contribution of the wall layer momentum (C < Cm) to the overall momentum 
balance are small. Under these assumptions we have : 

continuity equation 


dv dw 

+ dC 


= 0 


( 10 . 6 ) 
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momentum equation 


dv dv d ( p\ d ( t v \ 

V drj +W d( dri\Q) 3C\0/ 


The pressure is imposed by the external flow, so that 


(10.7) 


dp 

dC 


= 0 


The shear stress is determined from 


(10.9) 



(10.9) 


Using eqs(10.6) and (10.7) the integral form of the momentum equation is written 


^-/(c ! +pK = -- (10.10) 

an J Q 

' 0 


Neglecting the pressure term and substituting the velocity profiles from eqs(10.2) and 
(10.4) gives 


- -&c, 


( 10 . 11 ) 


In eq(10.10) the contribution to the integral for the region 0 < ( < Cm has been 
ignored since this is 0(Cm/&j), be., 0(1 /n) where n is large, particularly for large Reynolds 
number flows. 
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Equation(lO.ll) indicates that the momentum in the wall jet is reduced by the action 
of wall friction. However, for our purposes, the jet will travel only a very short distance 
before separating, so it is reasonable to assume that its momentum will remain constant. 


bjv 2 m = const (10.12) 

Having neglected curvature effects, eq(17) in ref.63 shows that the spreading rate of the 
jet will be constant, equal to that for a wail jet along a plane surface, i.e., 


where 


bj = Krj 


(10.13) 


K — 0.073 (10.14) 

is an experimentally determined constant. It follows from eq(10.12) that the velocity of 
the jet is reduced as 1/y/rf, i.e., 


— 



(10.15) 


where the blowing coefficient is defined as the ratio of the jet momentum to that of the 
external flield, just outside the boundary layer 


r - h nk 
M ~ Rvj 


(10.16) 
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10.4 Wall Jet Separation 


The only pressure gradient to which the jet is subject, after neglecting curvature 
effects, is the one due to the external flow. An approximate relationship for the influence 
of pressure gradient on wall shear stress for wall jets is 


- r „, 0 K'(m (^ dr ^j (10.17) 

where t^o is the shear stress at the wall with zero pressure gradient. From experiments by 
Bradshaw and Gee (ref.55) it is known that 


K' ~ 7 (10.18) 

4 

Thus separation (where r„ = 0) occurs when the following condition is satisfied 

teHIH (10 - 19) 

Table(lO.l) compares the separation criteria for the boundary layer and the wall jet. The 
right side is approximately the same for both cases. The wall jet, however, has greater 
momentum neax the wall. As a result, its characteristic dimension (distance of maximum 
velocity from the wall) is smaller than the corresponding characteristic dimension of the 
boundary layer (momentum thickness). In addition, higher velocities near the wall imply 
larger velocity gradients which result in greater shear stress at the wall. Thus, the first 
factor on the left side of the separation criterion is much smaller for the wall jet than for 
the boundary layer. As a consequence, the pressure gradient at separation is much larger 
for the wall jet and enables it to go farther against an adverse pressure gradient. 
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Using the following definitions 


Cf 0 = 


_ T v,0 




( 10 . 20 ) 


Cp 


P- Poo 

+ <X 2 ) 


( 10 . 21 ) 


t] = A d,R 


( 10 . 22 ) 


Re m 



(10.23) 


together with eqs(8.18b),(10.5),(10.13) and the experimental result (ref.63) 


Cf 0 = 0.0315i2e”° 182 ~ 0.004 ( 10.24) 

which is valid for Re — O(10 4 ), the separation condition eq( 10.19) transforms into 


A ^ 2 


21.527 [ V 2 

(i + « a ) kdc,/dn)\. 


c* 


(10.25) 


Equation(10.25) is plotted in fig(10.4). It is seen, that the blowing intensity required 
for a given displacement of the lower separation point depends only on the state of the 
boundary layer (i.e. whether it is laminar or turbulent), and is almost independent of the 
cone geometry and angle of attack, as is indicated by the almost horizontal curves. 
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10.5 Converged Solutions before and after Blowing 


Figures( 10.3a) and (10.3b) show the converged solutions for a cone with e = 5° at 
a = 30° for the two states of the boundary layer (laminar and turbulent) before and after 
blowing. The main observation, which reveals the beauty of the idea of blowing as a 
means of controlling separation, is that very small blowing intensities are required to move 
the separation points from their natural locations, as predicted by the viscous/ inviscid 
scheme in chapter 9, to points very close to the leeward generator. For both states of the 
boundary layer, blowing causes the separation to occur at a larger angle from the windward 
stagnation line, thus moving the vortices closer to the surface of the body toward the 
leeward generator. Smaller blowing intensity is required for the turbulent boundary layer 
for the same final configuration. This is explained by the fact that the separation for the 
turbulent boundary layer occurs naturally at a larger angle, and therefore the required 
A0, is smaller. 


10.6 Pressure Distribution 

The modified pressure distributions for the configurations shown in fig(10.3) are plot- 
ted in fig(10.5). It is seen, that blowing has the following effects : 

(i) It pushes the vortices (and as a result the vortex suction) closer to the leeward 
generator, thus closing the flow field. In the limit, as separation is suppressed completely, 
the results from the Jones theory are recovered. 

(ii) It reduces the size of the separation bubble. This is shown by the diminishing of 
the flat portion of the curves which represents the distance between the upper and lower 
separation points. 
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(iii) It weakens the vortices (as is shown from the diminishing vortex suction) by- 
pushing them closer to the surface. In effect, this reduces the vortex lift contribution, 
which is equivalent to reducing the angle of attack. Thus, it is seen that blowing allows 
control of the lift on a highly maneuverable aircraft without changing its attitude. 


10.7 Lift 

The relation between the lift and blowing coefficients is shown in fig(10.6). The fact 
that the curves drop more sharply as the relative incidence increases indicates that for a 
given body (e) blowing becomes more effective as angle of attack increases. 

The limit of all the curves is of course the Jones lift. Although it may seem impossible 
to eliminate separation completely, from fig(10.6) one may pick off each curve the point 
where the lift is within 5% of the Jones value. For all practical purposes, the vortex lift 
can then be neglected. 

As a reminder, it is repeated that at the high angles of attack to which some of the 
highly maneuverable aircraft operate, the main problem is to eliminate any asymmetries 
of the vortex system, vortex breakdown, or both. Thus, the desire to sacrifice some of the 
vortex lift in order to achieve this goal is not surprising. As Wood (ref.64) has pointed out, 
however, there is an exception to the rule that blowing reduces the vortex lift. This occurs 
when, for a given configuration, blowing stabilizes the vortex system which otherwise would 
have broken down. 
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11. EPILOGUE 


11.1 Discussion 

It is well known from previous studies that the “single line-vortex” model has the 
following disadvantages when used to represent the inviscid outer field about bodies at 
high angle of attack : 

(i) The position of the vortices is not very accurate. This should be expected, since 
the vortices are represented only globally in this model. To clarify this point a little further, 
we should remind ourselves that, in reality, the vorticity which is shed from the surface of 
conical bodies at incidence is distributed and not concentrated as the “single line- vortex” 
assumes. More complicated models which take this fact into account (see for example 
Smith ref.29) give vortex core locations which agree much better with experimental obser- 
vations. The predicted vortex locations are even worse when asymmetric vortex solutions 
are sought. Nevertheless, as Smith has pointed out (ref.8) the crude vortex locations given 
by the “single line-vortex” model are very useful as initial guesses for the more complex 
numerical “rolled-up core” model. 

(ii) The vortex lift is overestimated. This again is the result of a very strong suction 
generated on the upper surface of the body under the locations of the vortices. For delta 
wings, however, the non-linear lift is not a large part of the total unless the aspect ratio 
is very small, hence the error in the total lift is not too serious. The “rolled-up core” 
model also shows some suction, but the pressure peaks are much lower thus giving better 
agreement with experiments. 

(iii) Vortex solutions cannot be found below a minimum value of the relative incidence, 
which depends on the thickness of the wing and the location of separation. Experimental 
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observations (refs. 34-39), partially verify this result, since at small angles of attack the 
body radius, as it grows in the longitudinal direction, prevents the departure of free vor- 
tices. When the angle of attack becomes sufficiently high, the vorticity in the boundary 
layer accumulates along generators on the upper surface of the body. The vortices gener- 
ally do not separate from the body until some higher angle of attack is reached. However, 
the “rolled-up core” model gives solutions for much lower values of the relative incidence. 
The present work shows that, as long as the separation occurs on the upper surface, the 
minimum values of the relative incidence below which solutions do not exist for the single 
line-vortex” model are reasonable and agree well with experiments (see Bryson ref.20). 
Trouble occurs when solutions are sought for which the separation takes place on the lower 
surface of the body. First, the straight feeding sheet has to pierce through the wing in 
this case, and this is physically impossible. Second, as the angle of attack decreases the 
vortex approaches the separation point and inevitably comes around the leading edge in 
which case the feeding sheet assumes an almost horizontal position. When this happens, 
the force b alan ce between the vortex and its sheet is no longer possible, and as a result no 
solutions can be found. 

(iv) The pressure distribution is poorly predicted by this theory, principally because 
the vorticity in the feeding sheets is neglected. On the body surface, the pressure jumps at 
the point where the vortex sheet emanates. This is also physically impossible. In reality 
the vortex sheet adjusts its position and shape so that it coincides with a three-dimensional 
stream surface. Since the normal velocity across such a surface is zero, the force on the 
vortex sheet is zero as well. In this model, however, the pressure jump is necessary to 
create the force on the vortex sheet which balances the force on the vortex. 

Regarding the boundary layer solution for conical bodies the following may be said : 
the agreement of the predicted separation points with experiments is very good. Although 
this might have been expected when the terms dropped out of the cross-flow momentum 
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equation were found to be small, there was still the question of how an unrealistic pressure 
jump resulting from the “single line-vortex” model would affect the boundary layer solu- 
tion. Fortunately, because the lower boundary layer separates well before the point where 
the vortex sheet emanates (for the inviscid solution), the calculation of the boundary layer 
takes place in a region which is not affected much by the pressure jump across the vortex 
sheet. 


11.2 Conclusions 

(i) The “single line-vortex” model is limited in its accuracy but is adequate for the 
initial investigation of vortex flow control by tangential blowing. 

(ii) Displacement of the vortex separation has been shown to influence the location 
and strength of the vortices for both flat plate and elliptical cross-section conical bodies. 

(iii) The three-dimensional boundary layer over a circular cone has been analyzed. A 
method analogous to the Karman/Pohlhausen technique has been used to solve the cross- 
flow momentum equation, and the predicted separation lines agree well with experiments. 

(iv) Blowing tangentially from slots located symmetrically along cone generators near 
the point of cross-flow separation is an effective way to control vortex location and strength. 
For sufficiently large blowing the dependence on vortex lift can be drastically reduced, and 
the effects of flow asymmetries may be made negligible. 
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11.3 Recommendations for Further Research 


Additional work using an improved model should be undertaken in the following areas 
of the present model : 

(i) Inviscid outer field : A better representation of the vortex sheets is desirable 
in order to get more accurate vortex positions and eliminate the pressure jump on the 
surface of the body. In addition, inclusion of the secondary vortices which were mentioned 
in section(l.l) may indirectly affect the main vortex parameters by influencing the locations 
of separation due to their close proximity on the surface. 

(ii) Boundary layer : Two improvements are desirable in the boundary layer model. 
The first is an extension to non-circular cross-sections, and the second involves asymmet- 
rical vortex configurations. The first may be accomplished by approximating the metric 
coefficients for very thin elliptical cross-sections. For the second a two-parameter integral 
method is necessary in order to match the pressure at the edges of each separation bubble 
(right and left) simultaneously. 

(iii) Wall Jet : If blowing around leading edges with very small radius of curvature 
is desired (thin elliptical cross-section), then curvature effects must be included as is done 
in ref.63. Control of asymmetrical vortex shedding could also be analyzed in a similar 
manner provided that an appropriate boundary layer model is devised (see discussion in 
item (ii) above). 
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FIGURES 



WING 

WAKE 


Figure 1.1 Model of F-5F at a = 40°, in Northrop water tunnel. The vortex system 
is asymmetric, and the lower vortex has burst at some point over the wing. [G.E. Erickson, 
W.P.Gilbert : “Experimental Investigation of Forebody and Wing Leading Edge Vortex 
Interactions at High Angles of Attack” AGARD CP-342, No.ll, July 1983] 
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Figure 1.2 


Vortex formation over a slender delta wing at incidence. 



b. “Multiple line-vortex” model. 



c. “Single line-vortex” model. 



Figure 2.1 Three models representing vortex separation 
in the cross-plane of a conical body. 
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Figure 2.2 


Conformal mapping in the cross-plane. 




Figure 3.1 


A flat delta wing at incidence. 
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Figure 3.2 Schematic of the streamlines in the cross-plane 
of a flat delta wing at incidence with attached flow. 
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Figure 4.1 “Single-line vortex” model 

on a flat delta wing with leading-edge separation. 
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Figure 4.2 Locus of vortex positions 

for a flat delta wing with leading-edge separation. 
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Figure 4.3 Vortex strength vs. relative incidence 
for a flat delta wing with leading-edge separation. 
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Dimensionless Distance 



a. z-plane 



b. ("-plane 


Figure 4.5 


Contours of integration for the normal force. 
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Figure 4.6 Lift vs. relative incidence 

for a flat delta wing with leading-edge separation. 
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Figure 5.1 “Single-line vortex” model 

on a flat delta wing with displaced separation. 
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for a flat delta wing with displaced separation. 
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Figure 5.3 Vortex strength vs. relative incidence 
for a flat delta wing with displaced separation. 
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Figure 6.1 “Single-line vortex” model 

on a delta wing with elliptical cross-section. 



Figure 6.2 Schematic of the details in the cross-plane. 
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Figure 6.3 Loci of vortex positions 

for a delta wing with elliptical cross-section ( b/a = 0.05). 
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Figure 6.4 Vortex strength vs. relative incidence 
for a delta wing with elliptical cross-section ( b/a = 0.05). 
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Figure 6.8 Lift vs. relative incidence for delta wings with displaced separation. 



Fopple Curve 

(Locus of Vortex Positions for high cxj s) 



Figure 7.2 Loci of vortex positions on a circular cone. 
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«. typical streamline at outer edge of BL 

» typical streamline inside BL 

Figure 8.1 Secondary flow on a cone at incidence. 



(b) BL growth on a transversely curved surface. 
Figure 8.2 Streamline divergence producing a thinner BL. 
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vortices in separated region 




image vortices 

The mechanism of vortex coalescence within the BL. 



Figure 8.4 Coordinate system for the BL analysis on a circular cone 
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Figure 9.1 Schematic of BL separation 
in the cross-plane of a circular cone. 
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Figure 9.2 Flow chart for viscous/inviscid interaction. 
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Figure 9.3 Converged solutions for e = 5° and a - 30°. 
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; Figure 10.2 Wall jet profile. 
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BOUNDARY WALL 
LAYER JET 



Figure 10.3a Converged solutions before and after blowing 
(e = 5°, a = 30°, laminar BL). 
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Vortex 
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Figure 10.3b Converged solutions before and after blowing 
(s = 5°, a = 30°, turbulent BL). 
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Blown (C^ = 0.027, 0 S = 170°) 
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Figure 10.5b Pressure distribution before and after blowing 
(e = 5°, a — 30°, turbulent BL). 


Lift Coefficient 
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Figure 10.6 Lift vs. blowing. 


TABLES 
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Table 6.1 Summary of scaling laws for displaced separation. 


FLAT PLATE ELLIPSE (6/a = 0.1) 



vortex y-coordinate 


vortex z-coordinate 


vortex strength 


vortex lift 


yio (1 - fey /a) y 10 (l - fe y /a) 

Z 10 (l — fey/ a ) Zio 

k 0 (l - yjSy/a j k 0 (l - y/Sy/af" 00 

Clvo (l - \fSJa) 3 C LV 0 (1 - ^/<5 y /a) 2 - 6 


OVERALL LIFT IS OF THE FORM : 




£ \ £ 


m(b/a) 
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Table 8.1 Analogy between 2-D and conical BL. 


2-D 

CONICAL 

momentum 

equation 

v 2 (d0i/dx)+ 

(<5i + 29i)v e (dv e /dx) 

^(^22/^) + 

u e (^2 + 2 $ 22 )[(dv e /dr]) + (eu t /R)) 

first 
A = 

shape factor 

( 8 2 /v)(dv e /dx) 

( 8 2 /v)[(dv t /dri) + (euJR)] 

second 
K = 

shape factor 

ms 2 ) a 

m/s 2 ) a 

third 

H = }{K) 

shape factor 

{ 82 / 022 ) 

solution 

e 2 = 

(0.47i//v®) fV'dx 

{0A7v/ E L Vg) } E L v\dr} 

E = 

1 

exp[ 6 e](u e /v e )dri] 
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Table 10.1 Separation criteria for BL and wall jet. 
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APPENDICES 



APPENDIX 1 

Complex Potential for an Expanding Ellipse 


In ref.4, the complex potential for an expanding ellipse is given by 


/_n t , Uoo dS <7 + vV 2 -C* 

*'(*) ” u <* b °W + 2tt dx 2 


For an elliptical cone the following equations relate the geometrical variables 


S(x) = nab 


a = x tans ~ xs 


b = x tan 6 ~ x<5 


Differentiating eq(A1.2) and using eqs(A1.3) and (A1.4) yields 


dS 

— ~ 2nxs6 = 2naS = 2nbe 
dx 


The parameter bo(x) is defined by 
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(-41.1) 


(-41.2) 


(- 41 - 3 ) 


(41.4) 


( 41 . 5 ) 


where 


M*) = ao(x) In _ I / ^ 1„ (I - + 5 / ^ In (? ■ 


1 r da 0 


o+ 


<*e 


- ^a o (0 + )lnx - ia 0 (l )ln(l-x) 


.v 1 ds . 

° 0(X) = 2^^ 


dan 1 <£*5 


dx 2tt dx 2 


~ £(5 


For incompressible flow Moo = 0, so the first integral in eq(A1.6) becomes 


X 

ie«5 J ln(x — 0^ = -edx(\nx - 1) 


0 + 


The second integral in eq(A1.6) may be written as 


ie<5 J ln(f — x)d£ = -x)[ln(l -x) - 1] 


Finally, for the last two terms in eq(A1.6) we have 


a o (0 + ) = lim a 0 (x) = 0 


Thus, eq(Al.6) reduces to 


a 0 (l ) = lim a 0 (x) = eS 

r— ►!- 


bo(x) = —eS^x \n2yjx(l — x) — 1 + 2 } 


(Ai.s) 


Substituting now eqs(A1.5) and (A1.8) into eq(Al.l), the complex potential for an ellipse 
that expands in a conical manner is obtained 


= - Uoo e6 |x ln2^/a;(l - x) - 1 + ^}+ u °° &£ln 


<7 + \/cr 2 — c 2 


(A1.9) 


| 

I 

! 

I 

l 
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APPENDIX 2 

Distance around the Edge of an Ellipse 


The distance of the separation point from the leading edge of an elliptical cross-section 
may be expressed by the following integral 


77 = J dt] = J \J dy 2 + dz 2 


(A2.1) 


Here, y and 2 are related by 



(.42.2) 


, yb 2 

dz = -dy 

z a 2 


(A 2.3) 


2 2 

y a 


b 2 (a 2 - y 2 ) 


(.42.4) 


Substituting the above expressions into eq(A2.1) gives 


T) = 




b 2 



(A 2.5) 


Equation(A2.5) can be simplified, since a = 1, and the result is 
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(A2.6) 



i 


+ b 2 


l-y 


;dy 


or, since c 2 — a 2 — b 2 = 1 — b 2 , 


V = 



1 - c 2 y 

- y 2 


(A 2.7) 


In order to avoid the evaluation of the elliptic integral of the second kind, the integrand 
will be simplified further in the following manner 


1 ~ cV (1 4- cy)(l - cy) (l + c)(l -cy) (A2.8) 

l-y 2 (1 + y)(l - y) 2(1 - y) 

The last equality follows because close to the leading edge y — 1. 

Note that the above expression is valid for all thicknesses, i.e., from the flat plate case 
(c = 1) to the circle case (c = 0). Substituting eq(A2.8) into (A2.7) yields 


(A2 - 9) 

Integration by parts gives an approximate expression for the distance of the separation 
point from the leading edge of an elliptical cross-section, as a function of 8 y 


V = — 



(1 — C + C8y)8y + 


c — 1 
2 a/c 


In 


[C — C8y — 1 + ^/(l — C + C8y)c8y 
(1 - C)(l — C + C8y) 


(A2.10) 
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APPENDIX 3 

Evaluation of the Derivative dx/da 


For the term \s 2 in. eq(6.2) the derivative with respect to x can be evaluated directly from 

I 

I = ~ u ™ £ x { ln [ 2 V x(1 “ x) ] ” 2(1 - x) } 

For the rest of the complex potential ( Xcf + X*i)> we h ave 

i 

d\ _ dx dd dxcf dfh , dXc/ <$1 , dXcf dk dx c f dR dx,i db 

da ~ dd da dd\ da d6\ da dk da dR da db da 

\ 

and since 9\,Q\ (or < 7 a , <7i ), fc, a and b are all linear functions of x 

l 
1 

ddi _ 9i 
da a 


(A 3.2) 


(A3. 3) 


ddi 9\ 
da a 


(A3.4) 


dR _ R 
da a 


(A3. 5) 


db_b 

da a 
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(A3. 6) 


Equation(6.9) can be solved for the vortex strength 


k = 


(0, + 00(0, - 6 X ){9A ~ A 2 )(Mi + i? 2 ) 




\ / 


and differentiation with respect to a gives 


dk _ k 
da a 


(.43.8) 


In a similar manner, differentiating eqs(6.1) and (6.3) gives the following expressions 


dx 

50 


—lUooOt 




1 

0 + 01 + 


0i 

00 ! + R? 




(A3.9) 


?M- t k( 1 - - ^ + 1- 

ee, “1,0-0! 00i + i?v + 0i 


(A3. 10) 


dXcf 

50i 


= ik 




(A3. 11) 


d Xc f . /0i (0-0i)(00i+i? 2 ) \ 

5fc ^ ^ \0i (0 + 0i) (00i — R?') ) 


(A3. 12) 


dXcf _ 2RuooOi _ 0 ,p( 1 , i ^ fA3.13^ 


1 A 

oT = U«>£ln0 


(A3. 14) 


Finally, from fig(2.2) 


M _ d6_dc_ _ c 2 

da dc da 2a yj <r 2 — c 2 


(A3. 15) 


dd ^ M a ^ 
da 2 ^ \fa 2 — c 2 / 


(A3. 16) 
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APPENDIX 4 

Three-Dimensional Boundary Layer Equations 


A4.1 Equations 


Using an orthogonal system, which is not less general but simplifies matters consid- 
erably, the expression for a general element of length, is given by 


(*) 2 = h'ftdo 2 + h\(d n y + kifdo 2 


(-4.4.1) 


where the metric coefficients are in general functions of all three coordinates 


hi = M£,77,C) 
h 2 = h 2 (£, 77, C) 
h 3 = h 3 (£, t ], 0 

£ and 77 lie and are defined on the surface over which the boundary layer is flowing, while 
C extends into the layer. 

When the surface is regular, and not excessively curved in comparison with the bound- 
ary layer thickness 


hi = hi(Z,rj) 
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h 2 = h 2 (t,n) 


h 3 = h 3 (t,ri) 

If (£, rj, C) are known functions of some Cartesian system (x, y, z ), then 



(A4.2a) 


(A4.26) 


(A4.2c) 


The boundary layer equations in a general system of orthogonal curvilinear coordinates 
like the one just described, can be written as follows : 

continuity equation 


h 2 h 3 u ) + ^{h x h 3 v) + ^r( h x h 2 w ) = 0 


(A4.3) 


momentum equation in the £ - direction 


u v du w du uv dh\ 

h\ d£ + h 2 dr/ + h 3 d£ hih 2 drj 


v 2 

h\h 2 d£ 


1 dp 1 3rg 

3£ + 3C 


(A4.4) 


momentum equation in the 77 - direction 
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(-44.5) 


— ~ , v_dv w dv 
d£ h 2 dT) + h 3 d£ + 


uv dhj _ u 2 5/ij 

^2 ^ M 7 &T 


momentum equation in the C - direction 


Qh 2 dp gh 3 dQ 


dp 

37 = 0 

(44.6) 


For the derivation of the above equations, the flow has been assumed , , ■ 

with neither body forces nor Coriolis acceleration terms. stea dy, incompressible 

may also be written as Pressure gradient components 


J__dp 

Qh\ 


_Ue _ V^dUe U e V e dhj 

hl hdr) h 1 h 2 ~drf + 


vj dhj 
h\h 2 d£ 


(4 4.7) 


il dVg Vg dVg 

8h 2 5? ? *1 ^ ~ - 


!^dh 2 «2 aAi 

m 2 ae m 2 dp 


(44.8) • 


A4.2 Boundary Conditions 


(0 At the surface of the body (C - 0) th* « ,• „ 

y ^ the no-shp condition is 


u = v = w = 0 

(44.9) 


(' ) At the outer edge of the boundary layer (( — ► oo) 
the external flow 


the velocity should match that of 
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u -*• U e (t,T)) 


(A4.10<z) 


v-+v e (£,r]) (A4.106) 

w — * 0 (-4.4.10c) 


A4.3 Choice of the Coordinate System 

The equations in section (A4.1) are complicated mainly because of the presence of 
the metric coefficients hi,h, 2 ,fi 3 and their derivatives. It is therefore imperative that the 
coordinate system is chosen in such a way as to simplify both the differential equations as 
well as the boundary conditions. 

A first simplification applicable in boundary layer studies, is to restrict the general 
orthogonal system which was defined in section (A4.1), by setting 


h 3 U,T)) = 1 (44.11) 

which implies that ( represents an actual distance measured along a straight normal from 
the surface. As a result, only the choice of the two remaining surface coordinates £ and tj 
needs to be made. 

References 45-48 have an extensive discussion on the various possibilities for the choice 
of the two remaining coordinate axes. The problem which exists most of the times, is that 
there is usually one coordinate system in which the boundary layer equations take the 
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simplest form, and another one which offers the simplest boundary conditions. For the 
circular cone for example, since it is a developable surface (i.e., it can be rolled out into a 
plane without being stretched after suitable cuts have been made), a cartesian coordinate 
system exists such that h\ = h 2 = 1. However, none of the coordinates of this system 
lies along the cone generators so the simple “conical flow” boundary conditions are lost. 
Thus, the best choice seems to be an orthogonal coordinate system consisting of geodesics 
(surface curves connecting succesive points along the shortest route possible) and geodesic 
parallels. Then, the metric coefficient for the coordinates which are geodesics becomes 


/ll(£,7?) = 1 


(A4.12) 


and eq(A4.1) now reads 


(ds) 2 — (d£) 2 + h\(drj) 2 + (d() 2 (A4.13) 
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APPENDIX 5 

Solution of the Boundary Layer Equation 

A5.1 Laminar Boundary Layer 

For the laminar boundary layer the shear stress is given by 

dv dv 

If we define a dimensionless coordinate across the boundary layer 

A = - 

then it may be written 

= /(A) = aiA + 6jA 2 + ci A 3 -|- djA 4 for 0 < A < 1 

The constants Oi, 61, Cj, d\ will be evaluated from the boundary conditions. 

The first boundary condition eq(A4.9), applied at A = 0 reduces eq(A4.5) to 

dry, _ 1 dp 

d( dp 

Combining now eqs(A5.1),(A5.4) and (A4.8) yields 


(45.1) 


(45.2) 


(-45.3) 


(A5.4) 


- 142 - 


(A5.5) 


o 


d 2 v 


v 


dQ 


-!(Sr*" 


The first pressure gradient parameter may be defined in a manner analogous to that for 


the two-dimensional boundary layer 


A, 


P_ 

av 


dr] 


+ £U e 


(A5.6) 


Then eq(A5.5) reduces to 


d 2 v v 

u dC~~ Ve 6 2Av 


Evaluation of the second derivative on the left side by means of eq(A5.3) yields 


6, = - jA, (A5.7a) 

The second boundary condition eqs(A4.10), applied at A = 1 gives v = v e , dv/dC, — 
0 and d 2 v/d( 2 = 0. When these conditions are expressed in terms of eq(A5.3), yield 
repsectively 


/(l) = 1 =£• cti + + Cj + d\ — 1 


(A5.7b) 


/ , (1) = 0 ^ oq -+* 1b\ -(- 3ci -|- 4 d\ — 0 


(A5.7c) 
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/"( l) =0 =► b 1 +3c 1 + 6d 1 =0 


( Ab.ld ) 


Solving now the system of eqs(A5.7) gives 


_ 9 , 

a,-2+ — 


(A5.7e) 


Cl = 



(A5.7/) 


d i 



(A5.7g) 


Substitution of eqs(A5.7) into (A5.3) shows that 



(A5.8a) 


- = F(A) + A„G(A) 

V e 


(.45.86) 


where the functions F and G are exactly the same as in the two-dimensional flow 


F( A) = 1-(1- A) 3 (l + A) 


(A5.9a) 
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(A5.96) 


G(A) = |A(1 - A) 3 

Thus, the dimensionless quantities of interest will also have the same form as their two- 
dimensional counterparts 



(A5.10a) 


62 

6 


1 

= J(1-F- A„G)dA = 0.3 - 0.00833A„ 

0 


(A5.106) 


Oil 

8 



0.11746 


(A5.10c) 


0 - 22 

8 


1 

= J(F + A„G)(1 — F - A v G)d\ = 0.11746 - 0.00106A„ - O.OOOUA* 

0 


(A5.10d) 


812 

8 


1 

= J( 1 - F)(F + A„G)dA = 0.11746 - 0.0036A„ 
0 


(A5.10e) 


<5 



F - A v G)dX = 0.11746 - 0.0047A„ 


(A5.10/) 


Notice that A,, does not appear in the right side of functions which depend only on the 
axial flow, since there is no pressure gradient in the £ - direction. 
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The friction coefficient is now defined by 



(.45.11) 


from which it may be written 


TrjS _ T„ 6Vg 

fj.v e qv\ V 



(.45.12) 


Cross-flow' separation occurs when r,, = 0, or, when A, = —12, while stagnation points 
correspond to A,, = +7.052. 


Multiplying now eq(8.15) by {622! vv t ) yields 


v e&22 d9 2 2 , ^22 (ct 1 ^2 \ f & v e , __ ^ _ T r)fi @-22 
Rv di) + Rv V M [dr, 7 txv e 8 


(A5.13) 


which is identical in form with the corresponding equation for the two-dimensional bound- 
ary layer, eq( 10.26) in ref. 40. Following the same procedure for its solution as in the 
two-dimensional case, one may define the second pressure gradient parameter as 


r , 0L fdv e 

k ' l = ST ( W + £ “* 


Combining now eqs(A5.6) and (A5.14) gives 


K vL = 9 



( 45 . 14 ) 


(A5.14a) 
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which, when combined with eq(A5.10d) yields 


K vL = 


37 An ^ A, 


3969 l 5 15 144 , 


The boundary layer shape factor is defined by 


H 2 l = ^~ 


'22 


and using eqs(A5.10b)and (A5.10d) it can be written as 


(A5.146) 


(.45.15) 


H 2 l = 


0.3 - 0.00833A„ 


0.11746 - 0.00106A„ - 0.00011 A* 
Equations(A5.10d) and (A5.12) also combine to give 


= fi(K, L ) 


(A5.15a) 


t v 6 6 


22 


-(*♦£ 


[iv t 8 

Now eq(A5.13) can be written as 


— — — — — ^ — = f (K T 1 

5 15 144/ 63 ^ n ’ 


(A5.16) 


~ (H) °« + k ” l12 + = hk ' l) {A5 ' 17) 

and if the function F is introduced in a manner analogous to that for the two-dimensional 
case 


F(K vL ) = 2 f 2 (K vL ) - 4 K nL - 2K vL f 1 (K vL ) 


L( Z 1- 

63 V 5 15 


AT 

144, 


116 / 2 1 \ l2 
2 “ 315 A ’’ + V945 + 120/ A ” + 


- 1_ A 3 

9072 v 


(A5.18) 
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eq(A5.17) finally becomes 


d (ei 2 


F(K nL ) 


(A 5.19) 


dr] \Rv J v e 

which is a non-linear ordinary differential equation for (0 22 / Ru). The function F(K v l) can 
be approximated by a straight line 


F(K v l) = ci — diK v L 


(A5.20) 


where cl and di are constants. When this is done, eq(A5.19) transforms into 


dui 

dr) 


+ ui 


I, i \ 1 dv e u e 

(di — 1) — - — I- die — 
v e dr] v e 


= c L 


(.45.21) 


where 


U>L = 


av 


(^45.22j 


The integrating factor for eq(A5.21) is 


exp 


J PMdrj 

-0 


V 


d L -l 


El 


(A 5.23) 


where 


El = exp 


d L e j —dr] 
. J o Ve . 


(A 5.24) 
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and 




Then the solution to eq(A5.21) is written as 


uiL = 



o 


(.45.25) 


(A5.26) 


Using the dimensionless velocities defined in eqs(8.17) the solution may be expressed in 
the following way 


where 


6 2 


22 


RclU 1 
E lUoo V d <- 



(A5.27) 


El = exp 



(A5.28) 


The condition for separation in laminar flow is written in terms of the pressure gradient 
parameters 


A„ = -12 =► I< vL = -0.157 


(.45.29) 
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and from the definition of K v l in eq(A5.14) we get 


cl 

E l V*l 


' dV \ 11 

— + eUj j = -0.157 


(.45.30) 


The constants cl and di have the same values as for the two-dimensional laminar boundary 
layer (i.e., c L = 0.47 and d L = 6). These values may be substituted into eq(A5.30), and 
the separation criterion for laminar flow” can be expressed as 


SC L = El l V~ 6 



Ejy h dr] — - 0.334 


(A5.32) 


For sufficiently slender bodies and small angles of attack it may be assumed 


u e ~ lx,*, =*> U ~ 1.0 


(A5.33) 


so 


SC L = Ei l V~* 



-0.334 


(A5.32a) 


where 


El = exp 



(A5.2Sa) 
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To avoid infinite values of El at the reattachment points, the velocity will be approximated 
by a linear expression near these points 


V = ojql 7 ) where uj 0 l = const 


when this is done El becomes 


(A5.34) 


El = exp 



}dV_drj_ 
' V dV 


and since ( drj/dV = l/w 0 l = const the right side of the previous equation can be integrated 
to give 


£ l — y(6eM>t) 


(A5.35) 


Finally eq(A5.32a) gives 


SC L = \ (.A5.36) 

0 

which is the limit of SCl as rj approaches a reattachment point. 


A5.2 Turbulent Boundary Layer 


The procedure for the solution of the boundary layer equations for turbulent flow is 
similar to the one outlined in the previous section for the laminar case. The differences are 
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mainly due to the fact that the turbulent boundary layer grows thicker than the laminar 
one, as was explained in section (8.3) 

For the turbulent case it is necessary to define, as in the two-dimensional case, the 
following parameters 


Re v = 


v .6 


e< 7 22 


V 


(.45.37) 


k '^w. r4 {% +£U ’) 


(,45.38) 


Hi T = p- = 9 x(K,t) 
“22 


(A5.39) 


\c,„=^ = Rt-*9i(K, T ) 

2 QVj 


(A5.40) 


Substituting eqs(A5.37) through (A5.40) into (A5.13) yields after some algebra 


= F(K, T ) - (.45.41) 

where 


F(K v t) = 2.25^2 - (3.25 + 2.25 gi )I< vT (A5.42) 
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It must be noted that F{K v t) is again the same function as i 
and it may be approximated by a straight line 


F(K nT ) = c T - d T K v T 


When this is done, eq(A5.41) transforms into 


d.UT 

— (- u>t 

dr] 



u e , , 1 dv e 

e b d T —j- 

v e v e dr) 


where 



The integrating factor for eq(A5.44) is 


exp 


j Pt(i n)df) 




= E r vi T 


where 


Ej = exp 


£ ( dT + 1) / C 1 ’' 

Vtr 


and 


the two-dimensional case, 


(A5.43) 


Rct (A5.44) 


(A5.45) 


(A5.46) 


(A5.47) 
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(.45.48) 


P T (r,) = (d T + |) + 


dj dv e 
Ve dri 


Then the solution to eq(A5.44) is written as 


u>T — 


Rcj 


J E T vi T drj 


+ C2 


Wtr 


E T v d j 


(.45.49) 


or, using the same dimensionless terms defined by eqs(8.17) as well as c? = 0.016 and 
dj = 4 (same values as for the two-dimensional boundary layer) 


Et = exp 



(>15.50) 


9 


2.25 

22 


/ 1/ y- 25 0.016 R 
UoJ EtV 525 


■ rj 

[ J^rnV^/Jn 4- c. 

J ' 1 '* 

IHr 


( 45.51) 


The constant of integration c 2 may be determined by equating 0 2 2 from eq(A5.51) with its 
value for laminar flow, both being evaluated at the transition point. This yields 


c 2 = 26.729t/‘°- 125 £E 1 ' 125 ^ 9,5 



(.45.52) 


Separation of the turbulent boundary layer occurs when K v t = -0.06 as for the 
two-dimensionlal case. Using eqs(A5.38) and (A5.51) this condition translates to 
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= -3.75 


(A5.53) 


SC T = 


1 

E t V 5 


dV TT 

“7 h sU 

dr) 


m rj s 

J E T V 4 dri + c 2 

I)tr 


assuming again that U — 1, eqs(A5.47) and (A5.53) reduce to 


Et = exp 


5.25e J 

Vtr 


dr] 

V 


(4.5. 47a) 


SC T = 


1 

E t V 5 



■ Vs 

J ErV 4 dT} -t- c 2 

IHr 


-3.75 


(A5.53a) 


The linearity assumption for the velocity near reattachment points is employed again 


V = u 0 ti] where cjqt = const 


(A5.54) 


and the result is now 


£ t _ pr5.25(e/Wr) 


(-45.55) 


Substitution of the above approximations into eq(A5.53a) yields 


SCt = 


CJQX + £ 

5(o?ot + £) + 0.25s 


(A5.56) 


and since e is small, the second term in the denominator may be neglected leaving 
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(A5.57) 


scv = I 

which is the limit of SCt as 77 approaches a reattachment point for a slender cone, assuming 
that the boundary layer is turbulent from its start. 
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APPENDIX 6 
Program Listings 


This appendix contains five main FORTRAN programs : 

A. PROGRAM VORTEX-CIRCLE 

B. PROGRAM V-CP ELLIPSE 

C. PROGRAM V-CP CIRCLE 

D. PROGRAM K-CL 

E. PROGRAM BL 

The function of each program, as well as the function of each subroutine within the main 
progr ams , is explained with comments wherever is appropriate. 
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ooo ooooo 


PROGRAM VORTEX-CIRCLE 


C 

C THIS PROGRAM SOLVES THE SYSTEM OF EQUATIONS CONSISTING OF 

C THE SEPARATION AND FORCE-FREE CONDITIONS, FOR THE RIGHT 

C VORTEX POSITION. A CIRCULAR CROSS-SECTION IS CONSIDERED, 

C AND THE VORTEX SYSTEM CAN BE EITHER SYMMETRICAL OR 

C ASYMMETRICAL . 

C 

IMPLICIT REAL *8 (A-H,0-Z) 

DIMENSION H (10) , RHS (10) , WORK (10) ,X(10) , XOLD (10) , Z (10) 
DIMENSION RJAC (10,10) ,OA(10,10) ,IPVT(10) 

C 

C OUTPUT 

C 

c x IS THE SOLUTION OF THE SYSTEM IN VECTOR FORM 

C 

C RJAC 

C RHS 

C 

C H 

C ACC 

C INPUT 

OPEN (UNIT-6, FILE-' TTY : ' , STATUS-' NEW' ) 

1 CONTINUE 

NDIM IS THE DECLARED ROW DIMENSION OF THE JACOBIAN 

N IS THE ORDER OF THE MATRIX 

NMAX IS THE MAX ALLOWABLE # OF ITERATIONS 

N=4 

NDIM— 4 
NMAX— 25 
WRITE (5, 40) 

40 FORMAT ( ' GIVE THE RATIO ALFA/EPSILON ') 

READ <5, 100) AE 
WRITE (5, 41) 

41 FORMAT ( ' GIVE THE ASSUMED SEPARATION ANGLE ON THE RIGHT') 
READ (5, 100) SEPAR1 
WRITE (5, 42) 

42 FORMAT ( ' GIVE THE ASSUMED SEPARATION ANGLE ON THE LEFT') 
READ (5,100) SEPAR2 


INITIAL GUESS 


WRITE (5, 50) 

50 FORMAT ( ' GIVE Y10 ') 
READ (5,100) X(l) 

WRITE (5, 51) 

51 FORMAT ( ' GIVE Z10 ' ) 
READ (5, 100) X (2) 

WRITE (5, 52) 

52 FORMAT ( ' GIVE Y20 ') 
READ (5, 100) X (3) 

WRITE (5, 53) 


IS THE JACOBIAN OF THE SYSTEM 

IS THE RIGHT HAND SIDE OF THE EQ. : RJAC*H-RHS 
IN VECTOR FORM 

IS THE DIFFERENCE : XNEW-XOLD 

IS AN ESTIMATE OF THE MACHINE ACCURACY 
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53 FORMAT ( ' GIVE Z20 ') 

READ ( 5 , 100) X (4) 

100 FORMAT (F) 

C 

A=1 . 0D0 

GP=3 .141592654D0 

SEPl=SEPARl*GP/180 . 0D0 

SEP2- (180 . 0D0-SEPAR2) *GP/180 . 0D0 

C 

C ESTIMATE MACHINE ACCURACY 

C 

ACC=1 . 0D0 

140 ACC=0 . 5D0 *ACC 

ACCU-ACC+1 . 0D0 
IF (ACCU.GT . 1 . 0D0) GO TO 140 

C 

C PRINT INITIAL GUESS 

C 

DO 200 I-1,N 
WRITE (5,150) I,X(I) 

150 FORMAT ( ' X ( ' , II, ' ) , F10 . 3) 

200 CONTINUE 

C 

NITER-0 

NCONV-0 

C 

C LOOP FOR EACH ITERATION 

C 

250 CONTINUE 

DO 260 I— 1 , N 
XOLD ( I ) —X ( I ) 

260 CONTINUE 

C 

C FORM THE JACOBIAN MATRIX AND THE RHS VECTOR 

C 

CALL JACOB (NDIM,N,AE, SEP1 , SEP2 , X, RJAC, RHS ) 

300 FORMAT (IX, 4F15 . 8) 

C 

C SOLVE THE SYSTEM RJAC*H-RHS 

C 

CALL DECOMP (NDIM,N, RJAC, COND, IPVT, WORK, OA, Z) 

C 

C THE SYSTEM WILL BE SOLVED ONLY IF RJAC IS WELL CONDITIONED 

C 

COND I -COND +1 . 0D0 

IF ( COND I. EQ. COND) GO TO 390 

GO TO 420 

390 WRITE (5,400) 

400 FORMAT (' MATRIX IS SINGULAR TO WORKING PRECISION. ') 

GO TO 710 
420 CONTINUE 

CALL SOLVE (NDIM, N, RJAC, RHS, IPVT, H) 
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c 

c COMPUTE NEW VECTOR X 

C 

DO 450 1=1, N 
X ( I) =X (I) +H (I) 

450 CONTINUE 

NITER=NITER+1 

WRITE (6,470) NITER, (X(I),I=1,N) 

470 FORMAT (IX, 13, 8F9.5) 

C 

C ' CHECK FOR CONVERGENCE 

C 

DO 500 1=1, N 

IF (ABS ( (XOLD(I)-X(I) ) /XOLD(I) ) . LE . 1 . E-5 ) NCONV=NCONV+l 
500 CONTINUE 

IF (NCONV.EQ.N) GO TO 600 
IF (NITER. LT.NMAX) GO TO 250 
WRITE (5,520) 

520 FORMAT (' SOLUTION DOES NOT CONVERGE ') 

WRITE (5 r 521) 

521 FORMAT (' WITHIN THE SPECIFIED NUMBER OF ITERATIONS ' ) 

GO TO 710 

C 

C WRITE SOLUTION 

C 

600 WRITE (5,650) 

650 FORMAT (///,' THE SOLUTION OF THE SYSTEM IS : '//) 

DO 700 1=1, N 
WRITE (5,670) I,X(I) 

670 FORMAT ( ’ X ( ' , II, ' ) =' , F15 . 5) 

700 CONTINUE 

R1=DSQRT (X ( 1 ) **2+X (2) **2) 

R2=DSQRT (X (3) **2+X(4) **2) 

WRITE ( 5 , 7 0 5 ) R1 , R2 
705 FORMAT (2F15. 5) 

710 WRITE (5,750) 

750 FORMAT (/,' DO YOU WANT TO CONTINUE ? (1 FOR YES-0 FOR'NO)>',$) 

READ (6,*) I LOG 
IF ( ILOG . EQ . 1 ) GO TO 1 
STOP 
END 


SUBROUTINE SPLIT (AE, SEPl, SEP2, Yl, Zl, Y2, Z2,U,W, JF) 


C THIS SUBROUTINE SPLITS A COMPLEX FUNCTION INTO ITS 

C REAL AND IMAGINARY PARTS. 

C 

EXTERNAL F 

REAL*8 Yl, Y2, Zl, Z2, SEPl, SEP2,U, W 
COMP LEX *8 S1,S2,V,F 
S1=CMPLX (Yl, Zl) 

S2=CMPLX(Y2,Z2) 

V=F (AE, SEPl, SEP2 , SI, S2, JF) 

U=REAL (V) 

W=AIMAG (V) 

RETURN 

END' 
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SUBROUTINE JACOB (NDIM, N, AE, SEP1, SEP2, X, RJAC, RHS) 

C 

C THIS SUBROUTINE COMPUTES THE JACOBIAN OF THE GIVEN 

C SYSTEM OF EQUATIONS. 

C 

IMPLICIT REAL *8 (A-H,0-Z) 

DIMENSION RJAC (N,N) , RHS (N) ,X(N) 

Y1=X(1) 

Z1=X (2) 

Y2=X(3) 

Z2=X (4) 

DYZ=0 . 001D0 

CALL SPLIT (AE, SEP1, SEP2, Yl, Z1 , Y2, Z2, FR1, FI1, 1) 

CALL SPLIT (AE, SEP1, SEP2, Yl, Zl, Y2, Z2, FR2, FI2, 2) 

RHS<1)— FR1 

RHS (2)— FI1 

RHS (3) — FR2 

RHS (4)— FI2 

C 

C CALCULATE THE DERIVATIVES WITH RESPECT TO Yl 

C 

Y1N-Y1+DYZ 

CALL SPLIT ( AE , SEP 1 , SEP2 , YIN, Z 1 , Y2 , Z2 , FR1Y1 , FI 1Y1 , 1 ) 

CALL SPLIT (AE, SEP1, SEP2, YIN, Zl, Y2, Z2, FR2Y1, FI2Y1, 2) 

RJAC (1,1)” (FR1Y1-FR1) /DYZ 

RJAC (2,1)- (FI1Y1-FI1) /DYZ 

RJAC (3,1)- (FR2Y1-FR2) /DYZ 

RJAC (4,1)— (FI2Y1-FI2) /DYZ 

C 

C CALCULATE THE DERIVATIVES WITH RESPECT TO Zl 

C 

Z1N— Zl+DYZ 

CALL SPLIT (AE, SEP1, SEP2, Yl, Z1N, Y2, Z2, FR1Z1, FI1Z1, 1) 

CALL SPLIT (AE, SEP1, SEP2, Yl, Z1N, Y2, Z2, FP.2Z1, FT2Z1, 2) 

RJAC (1,2)— (FR1Z1-FR1) /DYZ 

RJAC (2,2)— (FI1Z1-FI1) /DYZ 

RJAC (3,2)— (FR2Z1-FR2) /DYZ 

RJAC (4,2) -(FI2Z1-FI2) /DYZ 

C 

C CALCULATE THE DERIVATIVES WITH RESPECT TO Y2 

C 

Y2N-Y2+DYZ 

CALL SPLIT (AE , SEP 1 , SEP2 , Yl , Z 1 , Y2N , Z2 , FR1Y2 , FI 1Y2 , 1 ) 

CALL SPLIT ( AE , SEP 1 , SEP2 , Yl , Z 1 , Y2N, Z2 , FR2Y2 , FI2 Y2 , 2 ) 

RJAC (1, 3) — (FR1Y2-FR1) /DYZ 

RJAC (2, 3) — (FI1Y2-FI1) /DYZ 

RJAC (3,3)— (FR2Y2 -FR2 ) /DYZ 

RJAC (4,3)— (FI2Y2-FI2) /DYZ 

C 

C CALCULATE THE DERIVATIVES WITH RESPECT TO Z2 

C 

Z2N-Z2+DYZ 

CALL SPLIT (AE, SEP1, SEP2, Yl, Zl, Y2, Z2N,FR1Z2,FI1Z2, 1) 
CALL SPLIT (AE, SEP1, SEP2, Yl, Zl, Y2, Z2N,FR2Z2,FI2Z2,2) 
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RJAC (1,4)= (FR1Z2-FR1) /DYZ 
RJAC (2,4)= (FI1Z2-FI1) /DYZ 
RJAC (3, 4) = (FR2Z2-FR2) /DYZ 
RJAC (4, 4) = (FI2Z2-FI2) /DYZ 


RETURN 

END 


COMPLEX FUNCTION F (AE, SEP1 , SEP2 , SI, S2, JF) 

C 

REAL* 8 SEP1, SEP2 

COMPLEX* 8 SI, S2 , CS1, CS2, SP1, SP2 , CSP1, CSP2, Ql, Q2 
COMPLEX* 8 A1,A2, A3, A4,B1,B2,B3,B4,C1,C2,C3,C4 
COMPLEX* 8 D1,D3,E1,E3,H1, H3,G1,G3,F1,F2 


SP 1=A*CMPLX (COS ( SEP 1 ) , S IN ( SEP 1 ) ) 
SP2=A*CMPLX (COS (SEP2) , SIN (SEP2) ) 

CS 1-CON JG (SI) 

CS2=CONJG(S2) 

CSP 1-CON JG(SPl) 

CSP2— CONJG (SP2) 

A1=1.0+A**2/SP1**2 

A2=1.0+A**2/SP2**2 

A3— 1 . 0+A**2/CSP2**2 

A4-1 . 0+A**2/CSPl**2 

Bl=l .0/ (SP1-S1) -CS1/ (SP1*CS1-A**2) 

B2=l .0/ (SP2-S1) -CS1/ (SP2*CS1-A**2) 

B3=l .0/ (CS2-CSP2) +S2/ (S2*CSP2-A**2) 

B4=l .0/ (CS2-CSP1) +S2/ (S2*CSP1-A**2) 

Cl-1 .0/ (SP1-S2) -CS2/ (SP1*CS2-A**2) 

C2=l .0/ (SP2-S2) -CS2/ (SP2*CS2-A**2) 

C3=l .0/ (CS1-CSP2) +S1/ (S1*CSP2-A**2) 

C4-1 .0/ (CS1-CSP1) +S1/ (S1*CSP1-A**2) 

Dl= (2 . 0*CS1-CSP1) /A 
D3- (SP2-2 . 0*S2) /A 
El=l . 0+A**2/Sl**2 
E3-1 . 0+A**2/CS2**2 
H1=-CS1/ (S1*CS1-A**2) 

H3=S2/ (S2*CS2-A**2) 

Gl-1 .0/ (S1-S2) -CS2/ (S1*CS2-A**2) 

G3-1 .0/ (CS1-CS2) +S1/ (S1*CS2-A**2) 

Ql- (A1*C2-A2*C1) / (B2*C1-B1*C2 ) 
Q2-(A3*C4-A4*C3) / (C3*B4-C4*B3) 

FI- (0, 1) * (Q1*H1-Q2*G1+E1) + (Dl-A/Sl) /AE 

F2-(0, 1) * (Q2*H3-Q1*G3+E3) + (D3+A/CS2) /AE 

IF (JF.EQ.l) F-Fl 

IF (JF.EQ.2) F=F2 

RETURN 

END 
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SUBROUTINE SOLVE (NDIM, N, A, B, IPVT,X) 

C 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A (NDIM, N) , B (N) , IPVT (N) , X (N) 

C 

C SOLUTION OF LINEAR SYSTEM, A*X=B 

C DO NOT USE IF DECOMP HAS DETECTED SINGULARITY 

C 

C INPUT . . 

C 

C NDIM - DECLARED ROW DIMENSION OF ARRAY CONTAINING A 

C 

C N - ORDER OF MATRIX. 

C 

C A = TRIANGULARIZED MATRIX OBTAINED FROM DECOMP 

C 

C B - RIGHT HAND SIDE VECTOR 

C 

C IPVT - PIVOT VECTOR OBTAINED FROM DECOMP 

C 

C OUTPUT . . 

C 

C X = SOLUTION VECTOR, X. 

C 

INTEGER KB, KM1, NM1, KP1, I, K, M 
REAL T 
DO 1 I-1,N 
X ( I ) =*B ( I ) 

1 CONTINUE 

C 

C FORWARD ELIMINATION 

C 

IF(N.EQ.l) GO TO 50 

NM1=N— 1 

DO 20 K-1,NM1 

KP1=K+1 

M-IPVT (K) 

T=X (M) 

X(M)«X(K) 

X(K)=T 

DO 10 I-KP1,N 
X(I)=X(I) +A(I,K) *T 
10 CONTINUE 

20 CONTINUE 

C 

C BACK SUBSTITUTION 

C 

DO 40 KB«1,NM1 

KM1=N-KB 

K=KM1+1 

X(K)=X(K) / A (K, K) 

T=- X(K) 

DO 30 I»1,KM1 
X ( I ) =X ( I ) +A ( I , K) *T 
30 CONTINUE 

40 CONTINUE 

50 X(1)=X(1)/A(1,1) 

RETURN 

END 
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c 


SUBROUTINE DECOMP (NDIM, N, A, COND, IPVT , WORK, OA, Z) 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


1 


c 


c 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A (NDIM, N) , IPVT (N) ,OA(NDIM,N) , WORK (N) , Z (N) 

DECOMPOSES A REAL MATRIX BY GAUSSIAN ELIMINATION 
AND ESTIMATES THE CONDITION OF THE MATRIX. 

USE SOLVE TO COMPUTE SOLUTIONS TO LINEAR SYSTEMS. 

INPUT . . 

NDIM * DECLARED ROW DIMENSION OF THE ARRAY CONTAINING A. 
N * ORDER OF THE MATRIX. 

A = MATRIX TO BE TRIANGULARIZED . 

OUTPUT . . 

A CONTAINS AN UPPER TRIANGULAR MATRIX U AND A PERMUTED 
VERSION OF A LOWER TRIANGULAR MATRIX I-L SO THAT 
(PERMUTATION MATRIX) *A = L*U 

• COND = AN ESTIMATE OF THE CONDITION OF A. 

FOR THE LINEAR SYSTEM A*X =B, CHANGES IN A AND B 
MAY CAUSE CHANGES COND TIMES AS LARGE 
IF COND+l.O .EQ. COND , A IS SINGULAR TO WORKING 
PRECISION. COND IS SET TO 1 . OE+32 IF EXACT 
SINGULARITY IS DETECTED. 

IPVT * THE PIVOT VECTOR. 

IPVT (K) * THE INDEX OF THE K-TH PIVOT ROW 
IPVT (N) = (-1) ** (NUMBER OF INTERCHANGES) 

WORK , Z . . THESE VECTORS MUST BE DECLARED AND 
INCLUDED IN THE CALL. THEIR INPUT CONTENTS ARE IGNORED. 

THEIR OUTPUT CONTENTS ARE USUALLY UNIMPORTANT. 

OA.. THE ORIGINAL N*N MATRIX 

THE DETERMINANT OF A CAN BE OBTAINED ON OUTPUT BY 
DET(A) = IPVT (N) * A ( 1, 1) * A(2,2) * ... * A (N, N) . 

REAL EK, T, ANORM, YNORM, ZNORM 
INTEGER NM1 , I , J, K, KP 1 , KB , KM1 , M 

DO 1 1*1, N 
DO 1 J-1,N 
OA(I, J)*A(I, J) 

CONTINUE 

IPVT (N) * 1 
IF (N.EQ.l) GO TO 80 
NM1 - N - 1 
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non oooooo ooo ooo oo 


COMPUTE 1-NORM OF A 

ANORM - 0.0 
DO 10 J=1,N 
T-0.0 

DO 5 1=1, N 
T=T+ABS (A (I, J) ) 

5 CONTINUE 

IF (T.GT. ANORM) ANORM-T 
10 CONTINUE 

GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 

DO 35 K=l, NM1 
KP1=K+1 


FIND PIVOT 
M=K 

DO 15 I=KP1,N 

IF (ABS (A (I, K) ) ,GT . ABS (A (M, K) ) ) M-I 

CONTINUE 

IPVT(K) = M 

IF (M.NE.K) IPVT(N) — IPVT (N) 

T = A (M, K) 

A (M, K) = A (K, K) 

A(K, K) = T 

SKIP STEP IF PIVOT IS ZERO 

IF (T.EQ.0.0) GO TO 35 

COMPUTE MULTIPLIERS 

DO 20 I=KP1 , N 

A ( I, K) = -A ( I , K) /T 

CONTINUE 


25 

30 

35 

C 

C 

C 

C 

C 

C 

C 

C 


INTERCHANGE AND ELIMINATE BY COLUMNS 

DO 30 J=KP1,N 
T=A(M, J) 

A (M, J) — A (K, J) 

A(K, J)=T 

IF (T.EQ.0.0) GO TO 30 

DO 25 I=KP1,N 

A (I, J)=A(I, J) +A(I,K) *T 

CONTINUE 

CONTINUE 

CONTINUE 

COND = (1-NORM OF A) * (AN ESTIMATE OF 1-NORM OF A- INVERSE) 
ESTIMATE OBTAINED BY ONE STEP OF INVERSE ITERATION FOR THE 
SMALL SINGULAR VECTOR. THIS INVOLVES SOLVING TWO SYSTEMS 
OF EQUATIONS, (A-TRANSPOSE) *Y = E AND A*Z = Y WHERE E 
IS A VECTOR OF +1 OR -1 CHOSEN TO CAUSE GROWTH IN Y. 
ESTIMATE = (1-NORM OF Z) / (1-NORM OF Y) 


- 165 - 


C SOLVE (A-TRANSPOSE) *Y - E 

C 

DO 50 K=1,N 
T=0 . 0 

IF (K.EQ.l)GO TO 45 
KM1-K-1 
DO 40 1=1, KM1 
T=T+A(I,K) *WORK (I) 

40 CONTINUE 

45 EK=1 . 0 

IF(T.LT.O.O) EK=-1 . 0 
IF(A(K,K) .EQ.0.0) GO TO 90 
WORK (K) =- (EK+T) /A (K, K) 

50 CONTINUE 

DO 60 KB=1,NM1 

K=N-KB 

T=0 . 0 

KP1=K+1 

DO 55 I=KP1,N 

T=T+A ( I, K) *WORK (K) 

55 CONTINUE 

WORK (K) =T 
M=IPVT(K) 

IF (M.EQ.K) GO TO 60 
T=WORK (M) 

WORK (M) “WORK (K) 

WORK (K) “T 
60 CONTINUE 

C 

YNORM-0 . 0 
DO 65 1=1, N 

YNORM=YNORM+ABS (WORK (I) ) 

65 CONTINUE 

C 

C SOLVE A*Z - Y 

C 

CALL SOLVE (NDIM,N, A, WORK, IPVT, Z) 
C 

ZNORM=0 . 0 
DO 70 1=1, N 
ZNORM=ZNORM+ABS (Z (I) ) 

70 CONTINUE 

C 

C ESTIMATE CONDITION 

C 

COND-ANORM* ZNORM/ YNORM 
IF (COND.LT.1.0) COND-l.O 
RETURN 
C 

C 1-BY-l 

C 

80 COND-l.O 

IF ( A ( 1, 1 ) . NE .0.0) RETURN 
C 

C EXACT SINGULARITY 

C 

90 COND=l . 0E+32 

RETURN 
END 
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PROGRAM V-CP ELLIPSE 


C 

C THIS PROGRAM COMPUTES THE VELOCITY AND PRESSURE DISTRIBUTIONS 

C ON THE SURFACE OF A CONICAL BODY OF ELLIPTICAL OR FLAT 

C CROSS-SECTION WITH SEPARATED FLOW. 

C 

COMPLEX CV1,CV2,CV3,CV4,CV5,CV6,G,G1,G2,Q,QA,Q1,Q2,QF, S, SP, SP1 
COMPLEX SI, S2, T, TA, Tl, T1A, T2 , T2A, TS, VCF, WT, WT1, WT2, WR, WB, WA, WQ 

C 

C DATA 

C 

OPEN (UNIT— 6, FILE-' TTY:' , STATUS-' NEW' ) 

WRITE (5,10) 

10 FORMAT ( ' GIVE THE CONE SEMIAPEX ANGLE ' ) 

READ (5,95) EPSILO 
WRITE (5,20) 

20 FORMAT ( ' GIVE THE THICKNESS RATIO ') 

READ (5,95) B 
WRITE (5,30) 

30 FORMAT ( ' GIVE THE ANGLE OF ATTACK ' ) 

READ (5,95) ALFA 
WRITE (5,40) 

40 FORMAT ( ' GIVE THE DISPLACEMENT OF THE SP ') 

READ (5,95) DY 
WRITE (5,50) 

50 FORMAT ( ' GIVE THE VORTEX POSITION ') 

READ (5,95) SI 
WRITE (5,60) 

60 FORMAT (/,' WHICH SURFACE ? (-1 FOR LOWER-+1 FOR UPPER) >',$) 

READ (6,*) SURF 

95 FORMAT (2F) 

A-1.0 

PI-3.14159 

C 

AE-ALFA/ EPSILO 
C-SQRT (A**2-B**2) 

R— ( A+B) /2 

ALF— ALFA*P 1/180.0 

EPS— EPSILO*PI/180 . 0 

C 

C WHEN THE SP IS ON THE UPPER SURFACE DZ > 0 

C WHEN THE SP IS ON THE LOWER SURFACE DZ < 0 

C 

DZ- (B/A) *SQRT (DY* (2*A-DY) ) 

X-A/EPS 

Y-0.0 

100 CONTINUE 

C 

C FOR THE UPPER SURFACE Z > 0 

C FOR THE LOWER SURFACE Z < 0 

C 

Z— SURF*B*SQRT (1- (Y/A) **2) 

S— Y+ (0, 1) *Z 
S2— CONJG (SI) 

C 

C FOR THE UPPER SURFACE G > 0 

C FOR THE LOWER SURFACE G < 0 

C 
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c 


c 


c 


150 


IF (B.EQ.0.0) SUR-SURF 
IF (B.GT.0.0) SUR-+1.0 
G=SUR*CSQRT (S**2-C**2) 

Gl-CSQRT (Sl**2-C**2 ) 

G2=CSQRT (S2**2-C**2) 

T= (S+G) / 2 
Tl= (Sl+Gl) / 2 
T2= (S2+G2) /2 

SP1=CSQRT ( (A-DY) **2+2*DZ* (A-DY) *(0,1) -DZ**2-C**2) 
SP=0 .5* ( A-DY+ (0,1) *DZ+SP1) 

Ql= (SP+T2) * (SP-T1). * (SP*T2-R**2) * (SP*T1+R**2) 
Q2=SP**2* (T1+T2) * (R**2-T1*T2) 

Q= (Q1/Q2) *ALF 

QF=1/ (T-Tl) -1/ (T+T2) +T1/ (T*T1+R**2) -T2 / (T*T2-R**2) 
TA=-C**2/ (2*A*G) 

T1A=T1 / A 

T2A=T2/A 

RA=R/ A 

BA=B / A 

QA=Q/A 

TS“ ( 1+S/G) /2 

WT=- (0,1) *ALF* (1+R**2/T**2) -(0,1) *Q*QF+B*EPS/T 
WT1- (0,1) *Q* (1/ (T-Tl) -T/ (T*T1+R**2) +1/T1) 

WT2- (0, 1) *Q* (1/ (T+T2) +T/ (T*T2-R**2) -1/T2) 


CVl-CLOG(T-Tl) 

CV1R=REAL(CV1) 

CV1I=AIMAG(CV1) 

IF (SURF. EQ. +1.0. AND. Y. GE. (1-DY) ) CV1I=*CV1I-2*PI 
DCV1I=ABS (CV1I-CV1IO) 

IF (SURF. EQ. +1.0. AND. DCV1I.GT. 4.0) CV1I=CV1I+2*PI 
IF (SURF. EQ. +1.0. AND. Y.GE. (1-DY) ) CV1I=CV1I-2*PI 
CV1-CMPLX (CV1R, CV1I) 

CV1IO-AIMAG (CV1) 

WRITE (5, 150) Y, CV1I 


CV2=CLOG (T+T2) 

CV3=CLOG (T*T1+R**2) 

CV4=CLOG(T*T2-R**2) 

CV5=CLOG(Tl) 

CV6=CLOG(T2) 

WQ— (0,1)* ( CV1 -CV2 +CV3 -CV 4 +CV 6 -CV5 ) 

WR=2* (0, 1) *R*ALF/T-2* (0,1) *Q*R* (1/ (T*T1+R**2) +1/ (T*T2-R**2) ) 
WB=EPS*CLOG(T) 

WA= WT * T A+ WT 1 * T 1 A+WT 2 * T 2 A+ WQ * Q A+ WR* RA+WB * BA 
VXS-EPS* (B/X) * (ALOG (2*SQRT (X* (1-X) ) ) -1/ (2* (1-X) ) ) 

VX=*EP S * REAL (WA) -VXS 

VCF-WT*TS 

VY=REAL (VCF) 

VZ— AIMAG (VCF) 

VEL-SQRT (VX**2+VY**2+VZ**2) 

CP= ALF * * 2 - 2 * VX- VY **2-VZ**2 
CPEE=CP/EPS**2 
WRITE (21, 150) Y,CPEE 
FORMAT (2F10 .3) 

Y=Y+0 .01 

IF (Y . LE . 1 . 0 ) GO TO 100 

STOP 

END 
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PROGRAM V-CP CIRCLE 


C 

C THIS PROGRAM COMPUTES THE VELOCITY AND PRESSURE DISTRIBUTIONS 

C ON THE SURFACE OF A CIRCULAR CONE WITH SEPARATED FLOW 

C FOR THE CASE OF AN ASSYMETRICAL VORTEX SYSTEM. 

C FROM THE OUTPUT, ONE CAN ALSO DETERMINE 

C THE POSITION OF THE STAGNATION, SEPARATION AND 

C RE ATTACHMENT POINTS . 

C 

COMPLEX Al, A2, A3, A4,B1,B2,B3,B4,C1,C2,C3,C4 
COMPLEX CS1 , CS2 , CSP 1 , CSP2 , Dl, D2 , D3 , D4 , D5 
COMPLEX DWA, DWS , DWS 1 , DWS2 , DWS 3,K1,K2,S,S1,S2,SP1,SP2 
COMPLEX WA, WA1 , WA2 , WA3 , WK1 , WK2 , WS 1 , WS2 , WCS 1 , WCS2 


DATA 


OPEN (UNIT-6,FILE=' TTY: ' , STATUS-' NEW' ) 

1 CONTINUE 

WRITE (5,10) 

10 FORMAT ( ' GIVE THE SEMIAPEX ANGLE ' ) 

READ (5,100) EPSILO 
WRITE (5,20) 

20 FORMAT ( ' GIVE THE ANGLE OF ATTACK ' ) 

READ (5,100) ALFA 
WRITE (5,30) 

30 FORMAT ( ' GIVE THE SEPARATION ANGLE ON THE RIGHT ' ) 
READ (5,100) SEPAR1 
WRITE (5,40) 

40 FORMAT ( ' GIVE THE SEPARATION ANGLE ON THE LEFT ' ) 
READ (5,100) SEPAR2 
WRITE (5,50) 

50 FORMAT ( ' GIVE THE RIGHT VORTEX POSITION ' ) 

READ (5,100) SI 
WRITE (5,60) 

60 FORMAT ( ' GIVE THE LEFT VORTEX POSITION ') 

READ (5,100) S2 
WRITE (5,70) 

70 FORMAT ( ' GIVE THE INTERVAL STEP ' ) 

READ (5, 100) DETA 
100 FORMAT (2F) 

PI-3.141592654 


EPS— EPS ILO*P 1/180 . 0 
ALF— ALFA*P 1/180.0 
SEPl-SEPARl*PI/180 . 0 
SEP2— (180 . 0-SEPAR2) *PI/180.0 
DET— DETA*PI / 1 8 0 . 0 


A-1.0 

SP1-A*CMPLX (COS (SEP1) , SIN (SEP1) ) 
SP2— A*CMPLX (COS (SEP2) ,SIN(SEP2) ) 
CSl-CONJG(Sl) 

CS2— CONJG ( S2 ) 

CSP1— CONJG (SP1) 

CSP2— CONJG (SP2) 


- 169 - 


ETA— 90.0 

150 ET=ETA*P 1/18 0.0 

Y«A*COS(ET) 

Z-A*SIN(ET) 

S=CMPLX(Y,Z) 

A1=1+A**2/SP1**2 

A2=1+A**2/SP2**2 

Bl=l/ (SP1-S1) -CS1/ (CS1*SP1-A**2) 

B2=l/ (SP2-S1) -CS1/ (CS1*SP2-A**2 ) 

Cl=l/ (SP1-S2) -CS2/ (CS2*SP1-A**2) 

C2=l/ (SP2-S2) -CS2/ (CS2*SP2-A**2) 

Kl-ALF* (A1*C2-A2*C1) / (B2*C1-B1*C2) 

K2-ALF* (A1*B2-A2*B1) / (C1*B2-C2*B1) 

WA2— 2* (0, 1) *K1*A/ (S*CS1-A**2) 

WA3-+2* (0, 1) *K2*A/ (S*CS2-A**2) 

WS1— +K1* (0,1) / (S-Sl) 

WS2=-K2*(0,1)/(S-S2) 

WCSl*+(0, 1) *K1*A**2/ (CS1* (S*CS1-A**2) ) 
WCS2— (0,1) *K2*A**2/ (CS2* (S*CS2-A**2) ) 

C 

D 1-CLOG (S) 

C 

D2-CLOG (S-S2) 

D2R— REAL (D2) 

021— AIMAG (D2) 

IF (ET.GE.SEP2) D2l-D2l+2*PI 
D2-CMPLX(D2R,D2I) 

C 

D3— CLOG (S-Sl) 

D3R— REAL (D3) 

D3I— AIMAG (D3) 

IF (ET.GE.SEP1) D3I=D3I+2*PI 
D 3— CMP LX (D3R, D3I) 

C 

D4-CLOG (S-A**2/CS2) 

D4R— REAL (D4 ) 

D4I— AIMAG (D4) 

DD4I-ABS (D4I-D4IO) 

IF (DD4I.GT.4.0) D4I-D4I+2*PI 
D4— CMPLX (D4R, D4I) 

D4IO— AIMAG (D4) 

C 

D5— CLOG (S-A**2/CS1) 

D5R— REAL (D5) 

05 I— AIMAG (D5 ) 

DD5I-ABS (D5I-D5IO) 

IF (DD5I.GT.4.0) D5I-D5I+2*PI 
D5-CMPLX (D5R, D5I) 

D5 10— AIMAG (D5) 

C 

WK1— (0,1) *(D3-D5) 

WK2— + (0, 1) * (D2-D4) 

C 

DWS2— Kl* (0, 1) * (1/ (S-Sl) -CS1/ (S*CS1-A**2) ) 
DWS3-+K2* (0, 1) * (1/ (S-S2) -CS2/ (S*CS2-A**2) ) 

C 


- 170 - 


WA1=2* (0,1) *ALF* (A/S) 

WA-WA1 +WA2 +WA3 +EPS*D1 
DWS1—ALF* (0,1)* (1+A**2/S**2) 

DWA-WA+WS1*S1+WS2*S2+WCS1*CS1+WCS2*CS2+WK1*K1+WK2*K2 
DWS-DWS1+DWS2+DWS3+EPS* (A/S) 

VX-EPS*REAL (DWA) 

VY-+REAL(DWS) 

VZ— AIMAG(DWS) 

CP»ALF**2-2*VX-VY**2-VZ**2 
CPEE-CP/EPS**2 
WRITE (6,200) ETA, CPEE 
200 FORMAT (2F10.3) 

ETA-ETA+DETA 

IF (ETA. LE. 90.0) GO TO 150 
WRITE (5,300) 

300 FORMAT (/,' DO YOU WANT TO CONTINUE ?(1 FOR YES-0 FOR NO)>',$) 

READ (6,*) I LOG 
IF (ILOG.EQ. 1) GO TO 1 
STOP 
END 
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PROGRAM K-CL 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


100 

200 


10 

20 


THIS PROGRAM COMPUTES THE VORTEX STRENGTH K 
AND THE LIFT COEFFICIENT CL 

FOR A CONICAL BODY WITH ELLIPTICAL OR CIRCULAR 
CROSS -SECTION r EXHIBITING LEADING EDGE SEPARATION. 


COMPLEX S , G1 , G2 , TH1 , TH2 , T , T1 , QK1 , QK2 , QK, FL 


DATA 


AE : relative incidence 

S : vortex position in the physical plane 

B : ellipse semi -minor axis 

(0 for a flat plate, 1 for a circle) 

PARAM : = separation angle (for a circle) 

= distance of separation from the leading edge 
(for an ellipse or a flat plate) 


OPEN (UNIT=6, FILE-' TTY: ' , STATUS-' NEW' ) 

A-l . 0 
P-3. 14159 
CONTINUE 
WRITE (5,200) 

FORMAT (/,' Enter AE, S (complex) , B, PARAM below.') 
READ (6,*) AE, S, B, PARAM 


JONES LIFT (linear with angle of attack) 


CL1— 2*P*AE 


VORTEX LIFT (non-linear with angle of attack) 


C-SQRT (A**2-B**2 ) 

R— (A+B) / 2 

Gl-CSQRT (S**2-C**2) 

G2-CSQRT (CONJG (S) **2-C**2) 

TH1-S+G1 

TH2— CON JG ( S ) +G2 

IF (B.LT.1.0) GO TO 10 

THS— PARAM 

THSR— THS*P/180 . 0 

DY-A* (1-COS (THSR) ) 

DZ=A*SIN (THSR) 

GO TO 20 
DY-PARAM 

DZ=B*SQRT (1- (A-DY) **2) 

Tl— A-DY+ (0,1) *DZ 

T— (Tl+CSQRT ( (A-DY) **2+2* (0,1) * (A-DY) *DZ-DZ**2-C*C) ) / 2 
QK1= (0 . 25*TH2**2+T*TH2-R**2) / ( (T+0 .5*TH2) * (0 . 5*T*TH2-R**2 ) ) 
QK2- (R**2+T*THl-0 . 25*TH1**2) / ( (T-0 .5*TH1) * (0 . 5*T*TH1+R**2) ) 
QK=T**2/ (T**2+R**2) * (QK1-QK2) 

VK-l/QK 

IF (B . EQ .1.0) GO TO 250 
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FL- (1+ ( A+B) / (A-B) ) *G1+ (1- (A+B) / (A-B) ) *S 
GO TO 260 
250 FL=2*G1 

260 CL2«(4*P*AE/ (QK*A**2) ) *REAL (FL) 

CL=CL1+CL2 
WRITE (5,300) VK, CL 
300 FORMAT (2F20.3) 

WRITE (5,400) 

400 FORMAT (/,' Do you want to continue ? (1 for yes - 0 for no) >' ,$) 
READ (6,*) I LOG 
IF (ILOG.EQ. 1) GO TO 100 
STOP 
END 
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PROGRAM BL 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


1 

3 


11 


12 


13 


14 


15 


16 


17 


18 


19 


20 


THIS PROGRAM INTEGRATES THE CROSS-FLOW BOUNDARY LAYER EQUATION 
FROM A STAGNATION (OR REATTACHMENT POINT) TO THE REGION WHERE 
THE FLOW IS EXPECTED TO SEPARATE, 

ON THE SURFACE OF A CIRCULAR CONE WITH SEPARATED FLOW, 

FOR EITHER SYMMETRICAL OR ASSYMETRICAL VORTEX SYSTEMS. 

IN ORDER TO SATISFY KELVIN'S THEOREM, THE VORTEX (K1-K2) 
SUGGESTED BY THE CIRCLE THEOREM FOR THE ASSYMETRICAL CASE 
HAS NOT BEEN INCLUDED. 

THE VELOCITY AND PRESSURE AT EACH POINT ARE 
CALCULATED BY SUBROUTINE VEL. 

FOR THE CASE OF A LAMINAR BOUNDARY LAYER 

THE SEPARATION POINT IS REACHED WHEN THE VALUE OF THE 

FUNCTION SC (SEPARATION CRITERION) IS -0.334 . 

FOR THE CASE OF A TURBULENT BOUNDARY LAYER 

THE SEPARATION POINT IS REACHED WHEN THE VALUE OF THE 

FUNCTION SC (SEPARATION CRITERION) IS -3.75 . 


COMPLEX S1,S2 
EXTERNAL El, SCI 


DATA 


OPEN (UNIT-6, FILE-' TTY:' , STATUS-' NEW' ) 
CONTINUE 
WRITE (5,3) 

FORMAT ( / 

READ (6, 

WRITE (5 
FORMAT ( 

READ (5, 

WRITE (5 
FORMAT ( 

READ (5, 

WRITE (5 
FORMAT ( 


' STATUS OF BL ? (0 FOR LAMINAR- 1 FOR TURBULENT) >', $) 

) NBL 

11 ) 

GIVE THE ANGLE OF ATTACK ') 

20) ALFA 

12) 

GIVE THE CONE SEMIAPEX ANGLE ' ) 

20) EPSILO 
13) 

GIVE THE SEPARATION ANGLE ON THE RIGHT ' ) 


READ (5,20) SEPAR1 


WRITE (5 
FORMAT ( 


WRITE (5 
FORMAT ( 


READ (5,20) SI 


WRITE (5 
FORMAT ( 


READ (5,20) S2 


WRITE (5 
FORMAT ( 


READ (5,20) ETA0 


WRITE (5 
FORMAT ( 


14) 

GIVE 


THE SEPARATION ANGLE ON THE LEFT ' ) 


READ (5,20) SEPAR2 


15) 

GIVE 


16) 

GIVE 


THE RIGHT VORTEX POSITION ' ) 


THE LEFT VORTEX POSITION ' ) 


17) 

GIVE 


STAGNATION POINT LOCATION ' ) 


18) 

GIVE 


THE STARTING POINT ' ) 


READ (5,20) ETA1 
WRITE (5,19) 

FORMAT ( ' GIVE THE 
READ (5,20) DETA 
PI-3.141592654 
FORMAT (2F) 


INTERVAL STEP ' ) 
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c 

ALF«ALFA*P 1/18 0.0 
EPS=EPSILO*PI/180 . 0 
SEP1— SEPAR1*PI/180 . 0 
SEP2= (180 . 0-SEPAR2) *PI/180 . 0 
DET=DETA*P 1/180.0 
ET1=ETA1*PI/ 180.0 
C 

C APPROXIMATION OF SC IN THE NEIGHBORHOOD OF 

C THE STAGNATION POINT 

C 

ETA-ETA0 

22 ETAN-ETA+DETA 
ETN-ETAN*P 1/180 . 0 
ET=ETA*P 1/180 .0 

IF (DETA.GT .0.0) GO TO 23 

IF (ETA.LT.ETA1) GO TO 25 

GO TO 24 

23 IF (ETA.GT.ETA1) GO TO 25 

24 E-1.0 

IF (NBL.EQ.0) SC-1.0/6.0 

IF (NBL.EQ.l) SC-1.0/5.0 

E0-E 

SC0-SC 

GO TO 26 


CALCULATE THE INTEGRAL FUNCTION : E 


SC1-SC0 

CALL ROME (ALF, EPS, SEP1, SEP2 , SI, S2 , ET1 , ET, EIN) 
IF (NBL.EQ.0) RC-6 
IF (NBL.EQ.l) RC-4.25 
E-EXP (RC*EPS*EIN) 


EVALUATE THE DERIVATIVE : dV/d(eta) 


CALL VEL (ALF, EPS, SEP1, SEP2, SI, S2,ET, V, CP) 
CALL VEL (ALF ,EPS,SEP1,SEP2,S1,S2, ETN, VI , CP 1 ) 
DVETA— (VI -V) /ABS (DET) 


EVALUATE THE DERIVATIVE : dCp/d(eta) 


DCP ETA— ( CP 1 -CP ) /ABS (DET) 


COMPUTE THE SEPARATION CRITERION 


CALL ROMSC (ALF , EPS, SEP1, SEP2, SI, S2, ET1, ET, SCIN, NBL) 
IF (NBL.EQ.0) M-6 
IF (NBL.EQ.l) M-5 
SC2— (DVETA+EPS) *SCIN/ (E*V**M) 

SC-SC1+SC2 

26 CALL VEL (ALF , EPS, SEP1, SEP2, SI, S2,ET,V, CP) 

C 

PARAM— 21 . 526*V**2/ ( (1+ALF**2) *DCPETA) 

C 

WRITE (5,30) ETA, SC, CP, PARAM 
30 FORMAT (4F15.3) 
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ETA-ETAN 

IF (DETA.GT .0.0) GO TO 35 
IF (ET.GT.SEP1) GO TO 22 
GO TO 38 

35 IF (ET.LT.SEP1) GO TO 22 

38 WRITE (5/ 40) 

40 FORMAT DO YOU WANT TO CONTINUE ? (1 FOR YES-0 FOR NO)>',$) 

READ (6,*) I LOG 
IF (ILOG.EQ.l) GO TO 1 
STOP 
END 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

101 


103 

c 

c 

c 


SUBROUTINE ROME (ALF, EPS , SEP1 , SEP2, SI, S2 , ET1, ET, RES) 


THIS PROGRAM COMPUTES THE INTEGRAL FOR THE FUNCTION "E" 
BY THE ROMBERG METHOD. 

THE INPUTS ARE : 

El : THE FUNCTION TO BE INTEGRATED 

ET1 : THE LOWER LIMIT 

ET : THE UPPER LIMIT 

ERR : THE DESIRED ACCURACY 

THE OUTPUT IS : 

RES : THE RESULT 


COMPLEX S1,S2 
EXTERNAL El 


"ZR” IS THE ARRAY OF APPROXIMATIONS 


DIMENSION ZR(10, 10) 


INITIALIZE THE INDEX AND COMPUTE THE FIRST APPROXIMATION 


1=1 

DEL-ET-ET1 

EI1-EI (ALF, EPS, SEP1, SEP2, SI, S2,ET1) 
EI2-EI (ALF, EPS, SEP1, SEP2, SI, S2, ET) 
ZR(1, 1) =0 .5*DEL* (EI1+EI2) 


THE MAIN LOOP. 

THE FIRST PART COMPUTES THE INTEGRAL USING A 2J+1 POINT 
TRAPEZOID RULE. THE METHOD MAKES MAXIMAL USE OF THE 
VALUES ALREADY COMPUTED. 


J=2**(I-1) 

DEL-DEL/ 2 
1 = 1+1 

ZR (I, 1) =0 . 5*ZR (1-1, 1) 

DO 103 K=1,J 
XR-ET1+ (2*K-1) *DEL 

EI3-EI (ALF, EPS, SEP1, SEP2 , SI, S2 , XR) 
ZR ( I, 1) =ZR (1,1) +DEL*EI3 
CONTINUE 


DO THE RICHARDSON EXTRAPOLATION 
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105 

C 

C 

C 


C 

C 

c 


110 

115 


DO 105 K-2, I 

ZR(I,K)-(4** (K-l) *ZR(I,K-1) 
CONTINUE 


-ZR(I-1,K-1) ) / (4** (K-l) -1) 


ERROR CONTROL 


ERR-0.001 

DIFF— ABS (ZR(I, I) -ZR(I, 1-1) ) 
IF (DIFF.LT.ERR) GO TO 115 


THE MAXIMUM NUMBER OF ITERATIONS ALLOWED IS 10 


IF (I.LT.10) GO TO 101 
WRITE (5,110) 

FORMAT (' MORE THAN 10 ITERATIONS REQUIRED, CHECK PARAMETERS') 
STOP 

RES-ABS (ZR (I, I) ) 

RETURN 

END 


SUBROUTINE VEL (ALF, EPS, SEP1 , SEP 2 , SI, S2, ET, V, CP) 


C THIS PROGRAM EVALUATES THE CROSS -FLOW VELOCITY COMPONENTS 

C ON THE CROSS-SECTION OF THE BODY 

C 

COMPLEX A1,A2,B1,B2,C1,C2 

COMPLEX CS1,CS2,CSP1,CSP2,D1,D2,D3,D4,D5 

COMPLEX DWA, DWS , DWS 1 , DWS 2 , DWS 3 , K1 , K2 , S , S 1 , S2 , SP 1 , SP2 

COMPLEX WA, WA1, WA2 , WA3, WK1, WK2 , WS1, WS2 , WCS1, WCS2 

w — ~ — ^ — — ^ ^ — — — — 

A-1.0 

PI-3.14159 

SP 1— A*CMPLX ( COS ( SEP 1 ) , S IN ( SEP 1 ) ) 

SP2— A*CMPLX (COS (SEP2) , SIN (SEP2) ) 

CSl-CONJG(Sl) 

CS2— CONJG ( S2 ) 

CSP 1-CON JG(SPl) 

CSP2-CON JG ( SP2 ) 

Y— A*COS (ET) 

Z— A*SIN (ET) 

S— CMPLX (Y, Z) 

Al— 1+A**2/SP1**2 
A2— 1+A**2/SP2**2 
Bl-1/ (SP1-S1) -CS1/ (CS1*SP1-A**2) 

B2-1/ (SP2-S1) -CS1/ (CS1*SP2-A**2) 

Cl-1/ (SP1-S2) -CS2/ (CS2*SP1-A**2) 

C2— 1/ (SP2-S2) -CS2/ (CS2*SP2-A**2) 

Kl— ALF* (A1*C2-A2*C1) / (B2*C1-B1*C2) 

K2-ALF* (A1*B2-A2*B1) / (B2*C1-B1*C2) 

WA1— +2* (0,1) *ALF* (A/S) 

WA2— 2* (0, 1) *K1*A/ (S*CS1-A**2) 

WA3-+2* (0, 1) *K2*A/ (S*CS2-A**2) 

WS1— +K1* (0,1)/ (S-Sl) 

WS2— K2*(0,1)/(S-S2) 

WCS1-+ (0, 1) *K1*A**2/ (CS1* (S*CS1-A**2) ) 

WCS2— (0, 1) *K2*A**2/ (CS2* (S*CS2-A**2) ) 
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c 

c 


c 


c 


c 


c 


D1=CL0G (S) 


D2=CL0G (S-S2 ) 

D2R=REAL(D2) 

D2I=AIMAG(D2) 

IF (ET.GE.SEP2) D2I=D2l+2*PI 
D2=CMPLX(D2R,D2l) 


D3=CL0G (S-Sl) 

D3R=REAL(D3) 

D3I=AIMAG(D3) 

IF (ET.GE.SEP1) D3I=D3I+2*PI 
D3=CMPLX(D3R,D3I) 


D4“CL0G ( S-A**2 /CS2 ) 
D4R=REAL(D4) 

D4I=AIMAG(D4) 

DD4I=ABS (D4I-D4IO) 

IF (DD4I.GT.4.0) D4I=D4I+2*PI 

D4“CMPLX(D4R,D4I) 

D4IO“AIMAG(D4) 


D 5 “CLOG ( S-A* *2 /CS 1 ) 
D5R-REAL(D5) 

D5I“AIMAG(D5) 

DD5I=ABS (D5I-D5IO) 

IF (DD5I.GT.4.0) D5I=D5I+2*PI 

D5=CMPLX(D5R,D5I) 

D5IO“AIMAG(D5) 


WKl=-(0, 1) * (D3-D5) 

WK2“+(0, 1) * (D2-D4) 

DWS1--ALF* (0,1) * (1+A**2/S**2) 

DWS2=- (0,1) *K1* ( 1/ (S-Sl)-CSl/ (S*CS1-A**2) ) 

DWS3=+ (0,1) *K2* (1/ (S-S2) -CS2/ (S*CS2-A**2) ) 

WA=WA1 +WA2 +WA3 +EP S *D 1 

DWA-WA+WS1*S1+WS2*S2+WCS1*CS1+WCS2*CS2+WK1*K1+WK2*K2 
DWS“DWS1+DWS2+DWS3+EPS* (A/S) 

VX-EPS*REAL (DWA) 

VY=+REAL (DWS) 

VZ— AIMAG(DWS) 

V=SQRT (VY**2+VZ**2) 

CP=ALF * * 2 -2 * VX- VY * *2 -VZ * * 2 

RETURN 

END 


FUNCTION El (ALF, EPS, SEP1, SEP2, SI, S2, ET) 


THIS PROGRAM CALCULATES THE INTEGRAND FUNCTION FOR ”E" 


COMPLEX S1,S2 

CALL VEL (ALF, EPS, SEP1, SEP2, SI, S2,ET, V,CP) 

EI-1.0/V 

RETURN 

END 
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SUBROUTINE ROMSC (ALF, EPS, SEP1, SEP2, SI, S2, ET1, ET, RES, NBL) 


THIS PROGRAM COMPUTES THE INTEGRAL FOR THE SEPARATION 
CRITERION BY THE ROMBERG METHOD. 

THE INPUTS ARE : 

SCI : THE FUNCTION TO BE INTEGRATED 
ET1 : THE LOWER LIMIT 

ET : THE UPPER LIMIT 

EPS : THE CONE SEMIAPEX ANGLE (IN RADIANS) 

ERR : THE DESIRED ACCURACY 
THE OUTPUT IS : 

RES : THE RESULT 


COMPLEX SI, S2 
EXTERNAL El, SCI 


"ZR" IS THE ARRAY OF APPROXIMATIONS 


DIMENSION ZR ( 10 , 10 ) 


INITIALIZE THE INDEX AND COMPUTE THE FIRST APPROXIMATION 


1=1 

DEL-ET-ET1 

SCI1-SCI (ALF, EPS, SEP1 , SEP2 , SI , S2, ET1, ET1, NBL) 
SCI2-SCI (ALF, EPS, SEP1, SEP 2, SI, S2, ET1, ET, NBL) 
ZR (1, 1) =0 . 5*DEL* (SCI1+SCI2) 


THE MAIN LOOP. 

THE FIRST PART COMPUTES THE INTEGRAL USING A 2J+1 POINT 
TRAPEZOID RULE. THE METHOD MAKES MAXIMAL USE OF THE 
VALUES ALREADY COMPUTED. 


j=2** ( 1-1) 

DEL-DEL/2 

1 = 1+1 

ZR (1, 1) —0 . 5*ZR (1-1, 1) 

DO 403 K— 1 , J 
XR»ET1+(2*K-1) *DEL 

SCI3=SCI (ALF, EPS, SEP1, SEP2, SI, S2, ET1, XR, NBL) 

ZR (1,1) — ZR( I, 1) +DEL*SCI3 

CONTINUE 


DO THE RICHARDSON EXTRAPOLATION 


DO 405 K— 2 , I 

ZR ( I, K) — ( 4** (K-l ) *ZR (I, K-l) -ZR (1-1, K-l) ) / (4**(K-1)-1) 
CONTINUE 


ERROR CONTROL 


ERR— 0 .001 

DIFF-ABS ( ZR (1,1) -ZR(I, 1-1) ) 
IF (DIFF.LT.ERR) GO TO 415 
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THE MAXIMUM NUMBER OF ITERATIONS ALLOWED IS 10 


IF (I.LT.10) GO TO 401 
WRITE (5,410) 

10 FORMAT (' MORE THAN 10 ITERATIONS REQUIRED, CHECK PARAMETERS') 

STOP 

15 RES=ABS ( ZR ( I , I ) ) 

RETURN 

END 


FUNCTION SCI (ALF, EPS , SEP 1 , SEP2, SI, S2 , ET1, ET, NBL) 


THIS PROGRAM COMPUTES THE INTEGRAND FUNCTION 
FOR THE SEPARATION CRITERION. 


COMPLEX S1,S2 
EXTERNAL El 

CALL VEL (ALF, EPS, SEP1, SEP2, SI, S2,ET, V,CP) 

CALL ROME (ALF, EPS, SEP1, SEP2, SI, S2,ET1, ET, EIN) 
IF (NBL.EQ.0) RC-6 
IF (NBL.EQ.l) RC-4.25 
E=EXP (RC*EPS*EIN) 

IF (NBL.EQ.0) M-5 
IF (NBL.EQ.l) M-4 
SCI=E*V**M 
RETURN 


